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Preface 


The  primary  goal  of  this  research  was  to  develop  an  accurate  model  to  analyze  an 
aperture  fed  stacked-patch  microstrip  antenna.  The  limited  bandwidth  performance  of 
microstrip  antennas  has  been  a  major  impediment  to  their  wider  application  and  new  de¬ 
signs  are  constantly  being  explored  to  solve  this  problem.  Both  aperture  fed  and  stacked 
patch  antennas  have  been  investigated  independently  and  shown  promising  improve¬ 
ments  in  bandwidth  performance.  I  hope  that  a  combination  of  these  designs  will  provide 
even  better  bandwidth  characteristics. 

I  chose  a  full-wave  analysis  based  on  the  mixed  potential  integral  equations 
(MPIE)  for  the  basis  of  this  model.  The  model  provides  a  complete  description  of  the 
near-fields,  including  surface  waves,  of  the  antenna.  The  Green’s  functions  are  calculat¬ 
ed  and  expressed  as  Sommerfeld  integrals.  Several  numerical  techniques  to  solve  the  in¬ 
tegrals  are  developed  and  tested.  All  algorithms  are  shown  to  provide  accurate  and  effi¬ 
cient  solutions.  When  used  with  the  methods  of  moments,  the  results  of  this  research 
should  be  useful  to  accurately  analyze  a  stacked-patch  antenna. 

I  would  like  to  thank  my  advisor,  Maj  Harry  Barksdale,  for  his  help  and 
enthusiasm.  It’s  always  easier  to  work  hard  on  a  project  when  others  show  an  active  in¬ 
terest  in  your  efforts.  I  would  also  like  the  thank  the  members  of  my  thesis  committee, 
Capt  Gregory  Warhola  and  Capt  Philip  Joseph  for  their  helpful  inputs.  I  am  extremely 
indebted  to  Eh".  Juan  Mosig  and  Dr.  Fred  Gardiol,  their  published  works  taught  me 
virtually  everything  I  know  about  microstrip  antennas.  I  must  also  thank  Dr.  David 
Pozar  for  providing  the  original  idea  from  which  this  research  grew.  Finally,  I  would 
like  to  thank  my  wife  Connie,  for  somehow  finding  it  in  her  heart  to  stay  married  to  me 
these  last  eighteen  months  when  I  spent  more  time  with  my  computer  than  with  her. 

James  B.  Nazar 
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Notation 


Ai2  Magnetic  vector  potential.  Boldface  defines  this  as  a  vector  quantity.  The  su¬ 
perscript  b  signifies  region  b.  The  number  Jl.  in  the  subscript  defines  this  as 
the  vector  potential  in  dielectric  Jb  for  a  source  on  interface ib. 

Ap  Aperture 

a^,  b^  X  and  y  dimensions  of  magnetic  charge  cells  on  the  aperture. 

^2’  ^2  ^  y  dimensions  of  electric  charge  cells  on  interface  2b  (patch  1). 

a^,  b^  X  and  y  dimensions  of  electric  charge  cells  on  interface  3b  (patch  2). 

a^,  bf  X  and  y  dimensions  of  electric  charge  cells  on  the  feedline, 

bjj^  Thickness  of  dielectric  la. 

bjjj  Thickness  of  dielectric  lb. 

^2b  Total  thickness  of  dielectrics  lb  and  2b. 

C Observer  cell  on  the  aperture,  i  is  an  index  variable. 

C21  Observer  cell  on  interface  2b  (patch  1). 

Observer  cell  on  interface  3b  (patch  2). 

Cj-j  Observer  cell  on  the  feedline. 

Nj  X  Nf  sub-matrix  used  to  calculate  in  the  i  cells  on  the  ground  plane 
due  to  electric  sources  in  the  y  cells  on  the  feedline. 

Nf  X  Nj  sub-matrix  used  to  calculate  in  the  i  cells  on  the  feedline  due  to 
magnetic  sources  in  the  j  cells  on  the  aperture. 

N|^  X  N2  sub-matrix  used  to  calculate  in  the  i  cells  on  the  ground  plane 
due  to  electric  sources  in  the  j  cells  on  interface  2b  (patch  1). 

[Cij'^]  Nj  X  Nj  sub-matrix  used  to  calculate  in  the  i  cells  on  the  ground  plane 
due  to  electric  sources  in  the  j  cells  on  interface  3b  (patch  2). 
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E 

EFIE 

F 

^bxx 

'JA32 


=b 

Ge21 


,-.bxx 

Gpil 


H 

HED 

HMD 

j 


N2  X  N  j  sub-matrix  used  to  calculate  in  the  /  cells  on  interface  2b  (patch 

1)  due  to  magnetic  sources  in  the  j  cells  on  the  aperture. 

X  sub-matrix  used  to  calculate  E^^*^  in  the  i  cells  on  interface  3b  (patch 

2)  due  to  magnetic  sources  in  the  j  cells  on  the  aperture. 

Electric  field  vector. 

Electric  Field  Integral  Equation 
Electric  vector  potential. 

Green’s  function  used  to  calculate  contribution  to  electric  field  from  electric 
surface  current  source.  Superscript  b  signifies  region  6,  xx  represents  contri¬ 
bution  to  j:-directed  (first  x)  electric  field  from  x-directed  (second  x)  current 
source.  The  number  32  in  the  subscript  defines  this  as  the  contribution  to  the 
electric  field  on  interface  ^b  due  to  a  source  on  interface  2.b. 

Green’s  function  used  to  calculate  contribution  to  electric  field  from  magnetic 
surface  current  source.  The  double  overbar  defines  this  quantity  as  a  dyadic. 

Green’s  function  used  to  calculate  contribution  to  magnetic  field  from  electric 
surface  current  source. 

Green’s  function  used  to  calculate  contribution  to  magnetic  field  from 
magnetic  surface  current  source. 


Green’s  function  used  to  calculate  contribution  to  electric  field  from  electric 
surface  charge. 


Green’s  function  used  to  calculate  contribution  to  magnetic  field  from  mag¬ 
netic  surface  charge. 

Magnetic  field  vector. 

Hankel  function  of  the  second  kind,  order  n. 

Horizontal  Electric  Dipole 
Horizontal  Magnetic  Dipole 
Imaginary  number  ^7J 
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Notation 


J2  Electric  surface  current  at  interface  2b. 

Electric  surface  current  at  interface  3b. 

Jf  Scattered  electric  current  on  the  feedline. 

Incident  electric  current  on  the  feedline. 

Jpi  Bessel  function  of  first  kind,  order  n. 

Complex  wave  number  for  dielectric  lb. 

M|  Magnetic  surface  current  on  the  ground  plane. 

Magnetic  surface  current  on  the  ground  plane  in  region  a. 

Magnetic  surface  current  on  the  ground  plane  in  region  b. 

MPIE  Mixed  Potential  Integral  Equation 
N  ^  Total  number  of  current  cells  the  aperture  is  divided  into. 

N2  Total  number  of  current  cells  patch  1  (interface  2b)  is  divided  into. 

Total  number  of  current  cells  patch  2  (interface  3b)  is  divided  into. 

NA  Not  Applicable 

Total  number  of  current  cells  the  feedline  is  divided  into. 

PI  Patch  1 

P2  Patch  2 

PV  Denotes  the  Cauchy  principle  value  of  an  integral  with  a  simple  pole, 

q  Electric  surface  charge, 

q^^  Magnetic  surface  charge. 

R  Radial  distance  between  source  and  observer  positions  (does  not  include  z- 

component  of  separation). 

S  jj  Source  current  cell  on  the  aperture,  /  is  an  index  variable. 
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Notation 


S2j  Source  current  cell  on  interface  2b  (patch  1). 

S^j  Source  current  cell  on  interface  3b  (patch  2). 

Sq  Source  current  cell  on  the  feedline. 

tjj  Thickness  of  dielectric  layer  2b  (b2jj  *  bjjj). 

T  Vector  rooftop  basis  function  associated  with  current  source  cell  S  j j 

V  Electric  scalar  potential. 

V  Electric  scalar  potential  for  an  electric  point  charge, 
s 

Magnetic  scalar  potential. 


^mq  Magnetic  scalar  potential  for  a  magnetic  point  charge. 
X,  y,  z  Unit  vectors  for  rectangular  coordinate  system. 


Nj  X  Nj  sub-matrix  used  to  calculate  in  the  i  cells  on  the  ground  plane 
due  to  magnetic  sources  in  the  j  cells  on  the  ground  plane. 

Nj  X  Nj’  sub-matrix  used  to  calculate  in  the  i  cells  on  the  feedline  due  to 
electric  sources  in  the  j  cells  on  the  feedline. 

^2  ^  ^2  sub-matrix  used  to  calculate  in  the  i  cells  on  interface  2b  (patch 
1)  due  to  electric  sources  in  the  j  cells  on  interface  2b. 

N2  X  N3  sub-matrix  used  to  calculate  E*^*^  in  the  i  cells  on  interface  2b  (patch 

1)  due  to  electric  sources  in  the  j  cells  on  interface  3b  (patch  2). 

X  N2  sub-matrix  used  to  calculate  E^^  in  the  i  cells  on  interface  3b  (patch 

2)  due  to  electric  sources  in  the  j  cells  on  interface  2b  (patch  1). 

N3  X  sub-matrix  used  to  calculate  E*^*^  in  the  i  cells  on  interface  3b  (patch 
2)  due  *0  electric  sources  in  the  j  cells  on  interface  3b  (patch  2). 
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Notation 


Column  vector  containing  the  amplitudes  of  the  electric  current  elements  on 
the  conductor  at  interface  2b  (patch  1). 


Column  vector  containing  the  amplitudes  of  the  electric  current  elements  on 
the  conductor  at  interface  3b  (patch  2). 


Column  vector  containing  the  amplitudes  of  the  reflected  electric  current 
elements  on  the  feedline. 


Column  vector  containing  the  amplitudes  of  the  incident  electric  current 
elements  on  the  feedline. 


[ctj  J  Column  vector  containing  the  amplitudes  of  the  equivalent  magnetic  current 
elements  on  the  aperture. 

V  Del  operator  in  transverse  coordinates  (;c,  y  or  p,  <p) 

Permittivity  of  dielectric  lb,  assumed  to  be  real  for  this  analysis. 

^bl2  Ratio  of  permittivitys  of  dielectric  lb  and  2b  (unitless). 

Critical  point  where  the  asymptotic  approximation  of  the  integrand  in  the 
Green’s  function  can  be  used. 

A.p  Location  of  pole  on  real  axis  for  integrands  in  the  Green’s  functions. 

rijj  Two  dimensional  pulse  doublet  basis  function  associated  with  electric  surface 
charges  in  source  cell  Sjj. 

p,  p'  Radial  position  vectors.  Unprimed  quantity  represents  observer  position  and 
primed  quantity  represents  source  position. 

p,  <|>,  2  Unit  vectors  for  cylindrical  coordinate  system. 

Permeability  of  dielectric  lb,  assumed  to  be  real  for  this  analysis. 

Pbi2  Ratio  of  permeabilitys  of  dielectric  lb  and  2b  (unitless). 
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Radial  frequency. 

Radial  cut-off  frequency  for  first  TE  surface  wave  mode  in  either  region. 

Angle  between  current  direction  and  direction  to  observer  point. 

Approximately  equal 
Multiply 


GREEN’S  FUNCTIONS  FOR  A  THEORETICAL  MODEL  OF  AN 
APERTURE  FED  STACKED-PATCH  MICROSTRIP  ANTENNA 


I.  Introduction 

The  concept  of  microstrip  antennas  was  first  proposed  by  Deschamps  back  in 
1953  [1].  However,  it  was  not  until  the  early  1970’s  that  practical  microstrip  antennas 
were  fabricated,  as  better  theoretical  models  and  photo-etch  techniques  for  copper  or 
gold-clad  dielectric  substrates  with  different  dielectric  constants,  low  loss  tangents,  and 
attractive  thermal  and  mechanical  properties  were  developed.  Microstrip  antennas  offer 
numerous  advantages  such  as  low  cost,  ease  of  construction,  thin  profile,  and  modular 
design.  They  can  easily  be  integrated  into  the  skin  of  an  aircraft  or  missile  without 
degrading  the  aerodynamics  of  the  vehicle.  Extensive  research  has  continued  up  to  the 
present  day  to  improve  microstrip  antennas  and  integrate  them  into  new  applications. 

1.1  Definition 

A  microstrip  antenna  is  composed  of  a  conducting  strip  radiator  separated  from  a 
ground  plane  by  a  dielectric  substrate.  The  input  to  the  antenna  patch  is  usually  supplied 
by  a  stripline  or  coaxial  probe  (see  Figure  1-1).  The  antenna  patch  conductors  arc 
normally  constructed  of  copper  or  gold  and  can  assume  virtually  any  shape,  but  are 
usually  rectangular  or  circular  shaped  to  simplify  analysis  and  performance  predictions. 
For  best  performance,  the  relative  dielectric  constant  of  the  substrate,  should  be  low 
(Ej.  =  2.5)  to  enhance  the  fringe  fields  which  account  for  the  radiation  [l]. 

1.1.1  Advantages  and  Disadvantages.  Microstrip  antennas  have  many 
advantages  over  conventional  antennas  and  can  be  used  in  several  applications  over  a 
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Antenna  Patch 


Plane 


Figure  1  - 1  Microstrip  antennas  with  (a)  coaxial  feed  and  (b)  microstrip  feed. 


frequency  range  of  approximately  100  MHz  to  50  GHz.  A  few  of  the  major  advantages 
of  microstrip  antennas  are: 

•  lightweight,  low  volume,  low  profile  planar  configurations  which  can  be 
made  conformal 

•  low  fabrication  cost;  readily  amenable  to  mass  production 

•  can  be  made  thin;  hence,  they  do  not  perturb  the  aerodynamics  of  host 
aerospace  vehicles 

•  the  antennas  may  be  easily  mounted  on  missiles,  rockets  and  satellites  without 
major  alterations 

•  the  antennas  have  low  scattering  cross  section 

•  linear  and  circular  (left  hand  or  right  hand)  polarizations  are  possible  with 
simple  changes  in  feed  position 

•  dual  frequency  antennas  easily  made 

•  no  cavity  backing  required 

•  microstrip  antennas  are  compatible  with  modular  designs  (Solid  state  devices 
such  as  oscillators,  amplifiers,  variable  attenuators,  switches,  modulators, 
mixers,  phase  shifters  etc.  can  be  added  directly  to  the  antenna  substrate 
board) 

•  feed  lines  and  matching  networks  are  fabricated  simultaneously  with  the 
antenna  [l] 
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Some  of  the  disadvantages  of  microstrip  antennas  are: 

•  narrow  bandwidth 

•  loss,  hence  somewhat  lower  gain 

•  most  microstrip  antennas  radiate  into  a  half  plane 

•  practical  limitations  on  the  maximum  gain  (=  20  dB) 

•  poor  endfire  radiation  performance 

•  poor  isolation  between  the  feed  and  the  radiating  elements 

•  possibility  of  exciting  surface  waves 

•  lower  power  handling  capability  [1] 

Narrow  bandwidth  severely  limits  the  applications  of  microstrip  antennas.  The 
bandwidths  of  the  microstrip  antennas  shown  in  Figure  1  -  1  are  typically  1-5%  of  the 
resonant  frequency.  Bandwidth  can  be  increased  by  increasing  the  thickness  of  the 
substrate  between  the  ground  plane  and  the  antenna;  but  loss,  radiation,  and  impedance 
mismatch  problems  arise  with  the  microstrip  or  coaxial  feeds  that  offset  any  bandwidth 
gains  [1]. 

1.2  Problem  Statement 

By  electromagnetically  coupling  a  microstrip  feed  on  a  separate  substrate  to  the 
microstrip  antenna  (Patch  1)  through  an  aperture  in  the  ground  plane  (see  Figure  1  -  2), 
the  substrate  thickness  can  be  increased  to  improve  the  bandwidth  while  avoiding  the 
above  mentioned  problems.  The  addition  of  a  parasitic  patch  (Patch  2)  overlaying  the 
antenna  patch  provides  additional  bandwidth  enhancement  as  well  as  additional  design 
control  over  the  antenna  radiation  pattern. 
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1 A  Research  Questions 

Answers  to  following  questions  will  be  found; 

1)  What  are  the  Green’s  functions  necessary  to  solve  for  the  various  currents  of 
the  antenna? 

2)  What  are  the  mathematical  characteristics  of  the  integrands  in  the  Green’s 
functions? 

3)  How  can  the  Green’s  functions  be  evaluated  numerically? 

4)  What  methods  can  be  used  to  reduce  computational  time? 

5)  How  can  the  solutions  of  the  Green’s  functions  be  used  in  a  moments 
method  solution  for  the  currents  of  the  antenna? 

1 J  Scope  and  Limitations 

To  make  the  exact  mathematical  analysis  of  the  microstrip  antenna  shown  in 
Figure  1  -  2  tractable,  certain  restrictions  are  made.  Each  dielectric  layer  is  isotropic, 
homogeneous,  and  lossless.  The  dielectric  layers  and  ground  plane  are  considered  infi¬ 
nite  sheets.  The  ground  plane  and  antenna  patch  conductors  are  infinitely  thin  and 
perfectly  conducting.  The  dielectric  layers  have  finite  thickness.  Expressions  are  found 
for  all  Green’s  functions  required  in  the  integral  equations  used  to  determine  the  electric 
and  equivalent  magnetic  currents  of  the  antenna.  The  Green’s  functions  consist  of  one  to 
three  integrals  that  cannot  be  evaluated  analytically.  Only  five  Green’s  functions  are 
explicitly  evaluated  to  demonstrate  the  various  numerical  integration  techniques  neces¬ 
sary  io  handle  the  different  characteristics  of  the  integrands.  These  techniques  can  then 
be  applied  to  all  other  Green’s  functions.  A  procedure  for  applying  a  moments  method 
solution  of  the  currents  of  the  antenna  is  only  described  in  this  paper,  no  actual  computa¬ 


tions  are  made. 
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1.6  Thesis  Organization 

The  remainder  of  this  document  is  organized  in  the  following  manner.  Chapter  II 
reviews  several  current  articles  on  aperture  fed  microstrip  antennas  and  methods  for  full- 
wave  analysis  of  microstrip  antennas.  Chapter  III  describes  the  solution  process  used  to 
obtain  the  Green’s  functions  and  the  moments  method  procedure  to  solve  for  the 
currents.  The  characteristics  of  the  integrands  in  the  Green’s  functions  are  explored  in 
Chapter  IV  and  the  necessary  numerical  integration  techniques  are  developed.  The  con¬ 
clusions  and  recommendations  of  this  research  are  discussed  in  Chapter  V.  There  are 
also  several  appendices  that  cover  background  information,  some  of  the  more  tedious 
calculations,  and  the  various  computer  programs  developed  during  this  research. 
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This  literature  review  outlines  current  research  on  microstrip  antennas.  Eight 
current  reports  on  aperture  fed  microstrip  antennas  and  the  means  used  to  analyze 
microstrip  antennas  will  be  summarized.  Since  an  exact  solution  of  the  currents  and 
electromagnetic  fields  of  the  proposed  antenna  is  most  desirable,  only  articles  dealing 
with  aperture  fed  microstrip  antennas  and  the  full-wave  analysis  of  microstrip  antennas 
will  be  reviewed. 


2 . 1  Aperture  Fed  Microstrip  Antennas 

D.  M.  Pozar  [2]  proposed  a  new  method  of  feeding  a  microstrip  antenna  where  a 
microstrip  antenna  on  one  substrate  is  fed  by  a  stripline  on  a  parallel  substrate  through  an 
aperture  in  the  ground  plane.  Pozar  notes  three  advantages  to  be  gained  from  an  apenure 
fed  configuration: 

(i)  The  configuration  is  well  suited  for  monolithic  phased  arrays,  where 
active  devices  can  be  integrated  on,  for  example,  a  gallium  arsenide 
substrate  with  the  feed  network,  and  the  radiating  elements  can  be  located 
on  an  adjacent  (low-dielectric  constant)  substrate,  and  coupled  to  the  feed 
network  through  apertures  in  the  ground  plane  separating  the  two 
substrates.  The  use  of  two  substrates  thus  avoids  the  deleterious  effect  of 
a  high-dielectric-constant  substrate  on  the  bandwidth  and  scan 
performance  of  a  printed  antenna  array. 

(ii)  No  radiation  from  the  feed  network  can  interfere  with  the  main 
radiation  pattern,  since  a  ground  plane  separates  the  two  mechanisms. 

(iii)  No  direct  connection  is  made  to  the  antenna  elements,  so  problems 
such  a  large  probe  self  reactances  or  wide  microstripline  (relative  to  patch 
size),  which  arc  critical  at  millimetre-wave  frequencies,  are  avoided.  [2] 


Pozar  used  a  simple  cavity  model  for  the  patch  antenna  and  small-hole  coupling  theory  to 
design  a  prototype.  Although  he  did  not  cite  specific  bandwidth  performance  results  for 
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the  prototype  antenna,  he  did  report  that  the  antenna  produced  a  normal  radiation  pattern. 

Sullivan  and  Schaubert  [3]  continued  with  Pozar’s  ideas  by  developing  an  exact 
mathematical  model  for  the  apenure  fed  microstrip  antenna  with  a  single  antenna  patch. 
They  developed  coupled  integral  equations  by  using  the  Green’s  functions  for  the 
grounded  dielectric  slabs  so  their  analysis  included  all  coupling  effects  and  the  radiation 
and  surface  waves  of  both  substrates.  The  analysis  was  eased  by  invoking  the 
equivalence  principle,  closing  off  the  aperture,  and  replacing  it  with  magnetic  surface 
currents  M^,  just  above  and  below  the  ground  plane  (see  Figure  2  -  1).  The  continuity  of 


y 


Figure  2-1  Antenna  and  feed  with  incident  and  induced  currents,  (a)  Original  problem,  (b) 

Equivalent  problem  [3). 


the  tangential  electric  field  through  the  aperture  was  maintained  by  making  the  magnetic 
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current  above  the  ground  plane  equal  to  the  negative  of  the  magnetic  current  below  it. 
With  the  space  below  the  ground  plane  (z  <  0)  denoted  as  region  a  and  the  space  above 
the  ground  plane  (z  >  0)  denoted  as  region  b,  the  electric  and  magnetic  fields  in  each 
region  were  written  as  a  summation  of  fields  due  to  the  various  currents: 

Er  =  Ea(Jinc)+Eajf)+Ea(Ms) 

Ha‘°'  =  Ha(Jinc}+Ha{Jf)+Ha{M,) 

Eb'°^  =  Eb(Jp)-Eb(M,) 

Hb'°^  =  Hb(Jp)-Hb(M,) 

where  the  known  incident  current  distribution  on  the  feedline  is  the  scattered 

current  on  the  feedline  is  Jf,  and  the  current  on  the  patch  is  Jp.  The  fields  on  the  right 
hand  side  of  the  above  equations  are  generated  by  the  specified  current  radiating  in  the 
presence  of  a  dielectric  slab  and  ground  plane  with  the  aperture  shorted. 

Coupled  integral  equations  are  obtained  for  the  three  unknown  currents  Jf,  Jp, 
and  Mg  by  enforcing  boundary  conditions.  A  Galerkin  moment  method  is  then  used  to 
solve  the  integral  equations.  The  solution  is  simplified  by  assuming  the  electric  currents 
on  the  antenna  patch  and  stripline  are  confined  to  the  y-  direction.  The  aperture  is 
assumed  electrically  short  and  the  magnetic  current  is  confined  to  the  x-  direction  (out  of 
the  page  in  Figure  2-1).  The  resulting  formula  for  Mg  is  not  quite  an  exact  solution 
because  the  formula  requires  one  parameter  to  be  determined  from  empirical  data.  The 
authors  compared  calculated  to  measured  results  for  input  impedance  for  several  different 
combinations  of  aperture  position,  dielectric  constant  and  dielectric  thickness  with  good 
agreement  [3]. 

The  most  recent  results  reported  in  the  literature  for  aperture  fed  microstrip 
antennas  are  from  Tsao  et  al  [4].  They  constructed  and  tested  an  aperture  fed  stacked- 
patch  microstrip  antenna  and  obtained  19.2%  bandwidth  for  input  VSWR  <  2  at  3.9  GHz. 
They  also  developed  a  two-input-port  feed  network  to  produce  either  dual  circular  or  dual 
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linear  polarization  modes,  again  with  approximately  20%  bandwidth.  Although  the 
authors  did  not  present  an  exact  analysis  for  their  antenna,  their  results  prove  that  an 
aperture  fed  stacked-patch  microstrip  antenna  can  yield  higher  bandwidth  than  a  more 
conventional  single  patch  microstrip  antenna. 

2.2  Analysis  of  Microstrip  Antennas 

To  find  the  electromagnetic  fields  radiated  by  a  microstrip  antenna,  the  integral 
equations  for  the  currents  on  the  feed  network  and  antenna  patch  must  be  solved  first. 
Mosig  describes  in  detail  the  mixed  potential  integral  equation  (MPIE)  as  applied  to 
microstrip  structures  that  could  be  very  useful  for  solving  these  currents  [5].  He  uses 
Green’s  functions  associated  with  the  scalar  and  vector  potentials  which  are  calculated  by 
using  stratified  media  theory  and  are  expressed  as  Sommerfeld  integrals.  The  MPEE  is 
numerically  stable  and  can  be  solved  with  efficient  algorithms.  The  MPIE  is  solved  in 
the  space  domain,  rather  than  the  spectral  domain,  to  help  keep  a  good  physical  insight  to 
the  problem.  The  author  presents  several  different  basis  and  test  functions  that  can  be 
used  for  the  method  of  moments  solution  of  the  integral  equations.  The  solution  rate  of 
convergence  for  the  different  basis  and  test  functions  are  compared  to  derive  the 
optimum  combination  for  fastest  convergence.  The  MPIE  includes  contributions  by  both 
surface  waves  and  radiation.  Multilayered  substrates  and  multiple  conductors  (stacked 
patches)  can  be  handled  by  making  suitable  modifications  of  the  Green’s  functions  and 
increasing  the  number  of  unknowns. 

An  earlier  paper  by  Mosig  and  Gardiol  provides  additional  insight  into  using  the 
MPIE  for  the  solution  of  the  microstrip  antenna  problem  [6].  In  this  paper,  the  authors  go 
into  the  precise  details  of  the  formulation  of  the  spatial  Green’s  functions  used  in  the 
MPIE.  They  also  derive  approximations  for  the  near-  and  far-field  solutions  for  the 
microstrip  antenna  and  point  out  the  significance  of  surface  wave  effects  in  the  solution. 
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Because  the  antenna  studied  in  this  thesis  will  consist  of  several  different 
dielectric  layers,  analytical  techniques  applicable  to  stratified  media  will  be  needed. 
Although  Mosig’s  paper  indicates  that  the  MPIE  can  be  used  for  stratified  media,  he 
derives  the  Green’s  functions  using  only  the  horizontal  electric  dipole  (HED)  and  point 
charge,  and  does  not  explicitly  derive  the  Green’s  functions  for  stratified  media  [5]. 
Kong  develops  integral  expressions  for  the  electric  and  magnetic  fields  for  both  the  HED 
and  horizontal  magnetic  dipole  (HMD)  in  both  infinite  and  semi-infinite  stratified  media 
[71.  The  Green’s  functions  for  stratified  media  for  both  the  HED  and  HMD  can  be  easily 
extracted  from  these  expressions.  The  Green’s  function  for  the  HMD  is  needed  to 
calculate  the  fields  radiated  by  the  magnetic  currents  used  to  close  off  the  aperture  in  the 
ground  plane  of  the  antenna. 

Nirod  and  Pozar  [8]  developed  a  method  to  calculate  the  two-dimensional  Green’s 
function  in  the  spectral  domain  for  a  current  element  between  any  two  layers  of  a 
multilayer  substrate.  Their  solution  is  obtained  by  solving  a  “standard”  form  containing 
the  current  element  between  any  two  layers  and  using  an  iterative  algorithm  on  the  “stan¬ 
dard”  form  to  find  the  solution  in  any  other  layer.  Their  analysis  draws  the  following 
conclusions; 

1)  The  numerical  solution  of  the  Green’s  function  for  a  point  different  from  the 
plane  of  the  current  element  exciting  it,  converges  much  faster  than  for  a  point  in  the 
same  plane  as  the  current  element. 

2)  The  numerical  solution  converges  faster  for  thicker  layers. 

3)  The  numerical  solution  for  a  structure  with  a  ground  plane  converges  more 
slowly  than  a  structure  without  a  ground  plane.  [8] 

Alexdpoulos  and  Jackson  report  the  effects  of  multiple  dielectric  layers  on 
microstrip  antennas  [9].  Using  a  superstrate  cover  over  a  microstrip  antenna  on  a 
grounded  substrate,  they  derive  the  integral  expressions  for  the  electromagnetic  fields. 
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They  then  show  how  substrate-superstrate  resonance  conditions  can  be  established  to 
maximize  antenna  gain,  radiation  resistance,  and  radiation  efficiency.  They  also  develop 
criteria  for  nearly  omnidirectional  H  and  E  plane  patterns.  The  relationships  of 
dielectric  constants,  dielectric  thicknesses,  and  antenna  placement  within  the  dielectric 
layers  on  gain,  radiation  resistance,  and  radiation  efficiency  are  presented  graphically  and 
in  great  detail. 

2.3  Summary 

Microstrip  antenna  research  continually  produces  more  versatile  and  useful 
designs.  The  work  by  Pozar,  Sullivan  and  Schaubert,  Tsao  et  al,  and  Alexopoulos  and 
Jackson  has  proven  that  aperture  fed  single  and  stacked-patch  microstrip  antennas  can 
produce  greater  bandwidth  performance  and  offer  additional  coiitrol  over  the  radiation 
patterns  compared  to  microstrip  or  coaxial  fed  antennas.  Analytical  methods  developed 
by  Sullivan  and  Schaubert,  Mosig  and  Gardiol,  Kong,  and  Nirod  and  Pozar  can  be  used 
to  calculate  rigorous  solutions  for  microstrip  antennas  in  stratified  media.  By  applying 
these  exacting  analytical  methods  to  aperture  fed  stacked-patch  microstrip  antennas,  it 
will  be  possible  to  define  frequency,  bandwidth,  and  radiation  characteristics  and  develop 
rules  and  criteria  for  the  efficient  design  of  these  antennas. 
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3.1  Overview  of  Analysis 

The  goal  of  this  analysis  is  to  find  the  induced  electric  currents  and  charge  distri¬ 
butions  on  the  antenna  patches  and  feedline  and  equivalent  magnetic  current  and  charge 
distribution  over  the  aperture  for  a  given  incident  current  on  the  feedline.  These 
distributions  can  then  be  used  to  determine  the  resonant  frequency,  input  impedance, 
bandwidth,  radiation  pattern  and  other  operating  characteristics  of  the  aperture  fed 
stacked-patch  microstrip  antenna.  The  various  currents  of  the  antenna  are  depicted  in 
Figure  3-1.  The  electric  currents  are  designated  as  for  the  incident  current  on  the 
feed  line;  is  the  scattered  current  on  the  feedline;  J2  is  the  current  on  the  patch  be¬ 
tween  the  frrst  and  second  dielectric  layer;  and  is  the  current  on  the  patch  between  the 
second  dielectric  layer  and  free  space.  By  using  the  equivalence  principle  [10]  the  aper¬ 
ture  can  be  closed  off  and  replaced  by  magnetic  surface  currents  Mj  just  above  and 
below  the  ground  plane.  Although  Figure  3-1  shows  the  electric  currents  in  the  x-direc- 
tion  and  magnetic  current  in  the  y-direction,  this  analysis  includes  all  currents  in  the  x-y 
plane  of  each  conductor  and  the  aperture.  Continuity  of  the  tangential  electric  field 
through  the  aperture  is  maintained  by  making  the  magnetic  current  above  the  ground 
equal  to  the  negative  of  the  magnetic  current  below  [3].  The  total  electric  and  magnetic 
fields  in  region  a  and  region  b  due  to  the  various  currents  are  given  by 


E.‘°‘  =  E.{Ji„c)  +  E.(Mi)  (3  -  1) 

H  =  H.  (Ji„c)  +  H.  (Jf)  +  H.  (M,)  (3  -  2) 

Eb'"*  =  Eb  (J2)  +  Eb  (J3)  -  Eb  (Ml)  (3  -  3) 

Hb‘°‘  =  Hb  (J2)  +  Hb  (J3) -  Hb  (Ml)  (3  -  4) 
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The  fields  of  the  problem  are  expressed  using  the  mixed  potential  integral  equa¬ 
tion  (MPIE)  and  are  solved  in  the  space  domain.  The  method  of  using  both  vector  and 
scalar  potentials  (mixed  potentials)  to  solve  scattering  and  antenna  problems  is  discussed 
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by  Harrington  [11]  and  applied  specifically  to  microstrip  antennas  by  Mosig  [5].  The 
electric  fields  due  to  the  electric  and  magnetic  sources,  respectively,  are  derived  from  the 
following  scalar  and  vector  potentials 

E  =  -jo)A  -  VV  (3 . 5) 

E  =  -lVxF 

e  (3-6) 

where  A  is  the  magnetic  vector  potential,  V  is  the  electric  scalar  potential,  F  is  the  elec¬ 
tric  vector  potential,  and  j  =  V^.  The  magnetic  fields  are  derived  from 

H  =lVx A 

V-  (3-7) 

H  =  -jtoF  -  VVm  (3  .  g) 

where  is  the  magnetic  scalar  potential.  Equations  (3  -  5)  and  (3  -  7)  assume  only 

electric  sources  are  present,  and  (3  -  6)  and  (3  -  8)  assume  only  magnetic  sources  are 

present.  The  vector  and  scalar  potentials  are  in  turn  expressed  using  the  corresponding 

Green’s  functions  as  superposition  integrals  of  the  charge  and  current  densities 


A(p)  =  |Ua(pIp')-  J{p')ds' 
V4p)  =  j^Gjplpjq(pjds' 
F(p)  =  |gf(pIp')-  M(p')ds' 


4P)  =  i  Gm(plp')qin(p')ds' 


(3-9) 


(3  - 10) 


(3-11) 


nn 


(3  -  12) 

where  the  dot  in  (3  -  9)  and  (3-11)  indicates  a  dot  product  of  the  dyadic  Green’s  func¬ 
tion  with  the  vector  surface  current.  The  vectors  p  and  p'  represent  the  radial  positions 
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of  the  observer  and  source  points,  respectively,  from  the  z-axis.  The  current  and  charge 


densities  are  related  through  the  continuity  equations 

Vj+jcoq  =  0  (3-13) 

VM+jcoqm=0  (3-14) 

Equations  (3  -  9)  through  (3  - 14)  are  used  in  (3  -  5)  to  (3  -  8)  to  describe  the  elec¬ 
tromagnetic  fields  in  the  different  regions  of  the  antenna.  Four  coupled  integral  equa¬ 
tions  are  then  obtained  for  the  four  unknown  currents  J^,  M^,  J2  and  by  enforcing  the 
boundary  conditions  of  Table  3-1.  The  other  boundary  conditions  of  the  antenna 
structure  (E^^  =  0  on  the  ground  plane,  conditions  on  the  normal  fields  at  the  dielectric 
interfaces,  etc.)  are  incorporated  in  the  construction  of  the  Green’s  functions. 


Table  3  •  1  MPIE  Boundary  Conditions 

1)  =  0  on  antenna  patch  1. 

2)  E^^  =  0  on  antenna  patch  2. 

3)  E^^  =  0  on  the  feedline. 

4)  is  continuous  through  the  aperture. 


The  first  step  in  developing  a  mathematical  model  for  the  antenna  is  to  calculate 
the  vector  potentials  from  which  the  various  Green’s  functions  are  constructed.  The 
MPEEs  satisfying  the  conditions  of  Table  3  - 1  arc  then  formulated  and  a  moments  meth¬ 
od  that  could  be  used  to  solve  the  current  and  charge  distributions  is  described.  Since  no 
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approximations  are  made,  this  solution  would  completely  describe  the  near-fields  of  the 
antenna.  A  time  dependence  of  is  assumed  throughout  this  analysis. 

3.2  The  Vector  and  Scalar  Potentials 

By  definition,  the  Green’s  functions  are  potentials  created  by  a  unit  source  which 
is  an  electric  or  Hertz  dipole  horizontally  located  within  one  of  the  microstrip  surfaces  or 
the  ground  plane.  With  a  linear  system,  the  superposition  principle  applies  and  the  po¬ 
tentials  of  any  finite  source  can  be  determined  by  representing  the  source  as  a  continuum 
of  elementary  dipoles  and  then  integrating  the  contributions  of  all  the  elementary  sources 
[12].  Mosig  and  Gardiol  thoroughly  develop  the  theory  and  method  of  Green’s  functions 
for  arbitrary  microstrip  structures  in  [6]  and  [12].  However,  their  Green’s  functions  are 
constructed  for  only  a  single  dielectric  layer  antenna,  with  only  electric  sources  on  the  di¬ 
electric-air  interface,  and  dielectric  permeability  fixed  at  the  free-space  value.  These  re¬ 
sults  are  extended  to  accommodate  the  needs  of  this  research  by  constructing  Green’s 
functions  for  a  structure  with  two  dielectric  layers  with  different  permittivities  and  per¬ 
meabilities;  with  electric  sources  at  the  dielectric-dielectric  and  dielectric- air  interfaces; 
and  with  equivalent  magnetic  sources  at  the  ground  plane-dielectric  interface. 

3.2.1  H ED  at  interface  2b.  The  analysis  of  this  structure  is  begun  by  determining 
the  fields  created  by  an  HED  along  the  jc-axis  at  the  interface  of  dielectric  lb  and  dielec¬ 
tric  2b  (see  Figure  3  -  2)  having  a  unit  moment  Idx  =  1  A  m.  The  resulting  Green’s 
functions  can  then  be  used  in  (3  -  9)  and  (3  -  10)  along  with  (3  -  5)  and  (3  -  7)  to  obtain 
the  E  and  H  fields  of  the  structure  for  an  arbitrary  distribution  of  sources  on  interface  2b. 

The  X  and  z  components  of  the  magnetic  vector  potential  in  each  dielectric  and 
free  space  can  be  expressed  in  the  form  of  Sommerfeld  integrals  [12] 
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a;i2(p) 


=  [h^Up) 


a,2  sinh(uibz)  dX 


4  Hf '(Xp) 


Az12(p)  =  cos  (j) I  Hi  ^Xp)  cosh(uibz)  dX 


a522(p)  =  j  Ho  ^Xp)[bx2  sinh(u2bz)  +  c^i  sinh(u2bz)]  dX 
A^2(p)  =  cos  <|>J  nf  \Xp)[b^  sinli(u2bz)  +  c^  sinh(u2bz)]  dX 


(3-15) 


(3-16) 


(3-17) 


(3-18) 
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A^2(p)  =  J  HS?\>-p)dx2  exp(-U3bz)  dX 
A^2(p)  =  cos  <^j  nf  \xp)  dz2  exp(-U3bz)  dX 


(3  -  19) 


(3  -  20) 


where  3x2.  a^,bx2.  bz2.  etc.  are  unknown  coefficients  to  be  solved  for,  p  is  the  radial 
distance  Ipl  between  the  source  HED  and  the  observer  point  and 

uib  =  V  X^  -  k?b .  U2b  =  V  X^  -  k^b .  U3b  =  V  X^  -  kfb 

with  k^jj,  k2{j,  and  k^jj  being  the  wave  numbers  in  dielectrics  lb,  2b,  and  3b,  respectively 
(k^jj  is  assumed  to  be  free  space,  k^,  but  is  referred  to  here  as  k^jj  for  consistency  of 
notation).  The  wave  numbers  are  defined  as 

kib  =  coVeibPib.  k2b  =  o)Ve2b  P2b.  ksb  =  coVcb  P3b 
where  eib2b3b  ^^Ib2b3b  permittivity s  and  permeabilitys  of  the 

corresponding  medium.  Axi2  is  defined  as  the  x-directed  magnetic  vector  potential  in 
region  b,  dielectric  lb  (0  ^  z  <  b^jj)  for  a  HED  at  interface  2b.  The  other  values  of 
Ax,z  nm  are  defined  similarly.  A  definition  of  the  integration  path  C  along  with  some 
background  behind  the  development  of  (3  -  15)  -  (3  -  20)  can  be  found  in  Appendix  A. 

The  boundary  conditions  to  solve  for  the  coefficients  of  (3-15)  -  (3  -  20)  are 
developed  in  Appendix  B.  The  pertinent  results  are  summarized  as  follows: 

@z  =  bib 


Ax12  -  A^2 

(3-21) 

1  8A^12  1  9A^22 

L-«(P) 

Plb  P2b 

2np 

(3  -  22) 

Ab  Ab 

^z\2  _  Az22 

Plb  P2b 

(3  -  23) 
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(3)  z  =  b2b 


1 

1 

3Ay22 

ElbHlb 

dz 

e2bl^2b 

dz 

Ax22  = 

^^2 

1  ^A^22  _ 

1  3Ax32 

P2b  5z 

P3b  dz 

_ 

A^32 

P2b 

P3b 

1 

3A^2 

-  1 

^A^2 

e2bP2b 

3z 

e3bP3b 

dz 

/  1  1 

ie2bH2b  eib^ibl  5p 


^::r^  =  cos  <j)  ( - ^ 

dz  |e3bP 


1  I 

I^Ax32 

E2bP2bJ 

1  dp 

and  the  Sommerfeld  radiation  condition 


lim  r(~  +  jkA)  =  0 

r  — >oo  ^  ^  ' 


where  A  also  satisfies  the  homogeneous  Helmholtz  equation 

(v%k^)A  =  0 

The  Sommerfeld  condition  determines  the  branch  interpretation  [9]  of  U3b 

Vx^-k^b  1.  W^k3b 


U3b  = 
U3b=j 


3b 


U1  ^  k3b 


The  term  5(p)/(2jtp)  in  (3  -  22)  represents  the  unit  current  source  {Idx  =  1  y 
interface  lb,  where  5(p)  is  the  Dirac  delta  function  and  can  be  expressed  as  [12] 


(3  -  24) 

(3  -  25) 

(3  -  26) 

(3  -  27) 

(3  -  28) 

(3  -  29) 

(3  -  30) 
-  k3b  as 

\  m)  at 

(3-31) 
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Using  the  boundary  conditions  defined  above  and  (3  -  15)  -  (3  -  20),  eight  equations  are 
formed  and  solved  for  the  eight  unknown  coefficients  to  yield 


' 

Axi2{p)  =  ^  h{)^\^p)  ^  sinh(uibz)  dX 

J  C 


(3  -  32) 


A2i2(p)  =  ^  cosh{uibz)  dX 


4k 


(3  -  33) 


A"x22(P)  =  ^| 


Nbx2(^)  sinh(u2bz)  +  Ncx2(^)  sinh(u2bz) 
De(x) 


dX 


(3  -  34) 


‘‘  HfW) 


A”r22(p)  = 


Nbz2(^)  sinh(u2bz)  -f  Ncz2(^)  sinh(u2bz) 


dX 


(3  -  35) 


A^x32(P)  =  ^(  exp(-u3bz)  dX 


(3  -  36) 


=  H^xp) 

where  N  represents  the  numerators  of  3*2,  a^,bx2,  b^,  etc.  and  D  represents  the 
denominators.  The  zeros  of 

D^X)  =  +  uib  coth(bibUib)]  U2b  cosh(u2b(b2b  *  bib)) 

'  +  [pbi2U^b  +  Pb23UibU3b  coth(bibUib)]  sinh(u2b(b2b  -  bib))  (3  -  38) 


exp(-u3bz)  dX 


(3-37) 
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\^[ebi3U3b  +  Uib  tanh(bibUii,)]  U2b  cosh{u2b(b2b  -  bib)) 

+[ebi2uib  +  eb23UibU3btanh(bibUib)]sinh(u2b(b2b-bib))  (3-39) 

define  the  surface  wave  poles  in  the  dielectrics  where 


and 


Ebll  ■ 


|J-bl2  = 


-Ilk 

e2b 

eb23  =  1^. 
£3b 

£bl3=f^ 

£3b 

M-lb 

II  -^2b 

Pb23  -  — , 
M-3b 

II  ^1*’ 

Pbl3  =  „ 

P-3b 

P2b’ 

The  roots  of  and  represent  the  frequencies  at  which  TE  and  TM  surface 
waves  propagate  [9].  The  terms  and  are  related  to  the  reflection  factors  on  a 
microstrip  structure  for  an  incident  perpendicular  polarized  wave  (TE)  and  parallel 
polarized  wave  (TM),  respectively  [12].  Surface  waves  represent  power  that  is  propagat¬ 
ed  along  the  surfaces  of  the  dielectric  layers,  rather  than  radiated  into  space.  The 
numerator  terms  of  (3  -  32)  -  (3  -  37)  are  given  in  Appendix  C. 

3.22  HED  at  interface  3b.  In  this  case  the  HED  is  at  the  interface  of  dielectric 
2b  and  dielectric  3b  (free  space)  as  shown  in  Figure  3  -  3. 

The  X  and  z  components  of  the  magnetic  vector  potential  in  each  dielectric  and 
free  space  can  be  expressed  as 


ax3  sinh(uibz)  dX 


Az13(p)  =  cos  (p 


»|hHXp) 


a^  cosh(uibz)dX 


A^3{p)  =  j  Ho  \Xp)  [b$3  sinh(u2bz)  +  0^3  sinh(u2bz)]  dX 
a523(p)  =  cos  <pj  Hf  ^Xp)[b|3  sinh(u2bz)  +  c^  sinh(u2bz)]  dX 


(3-40) 


(3-41) 


(3-42) 


(3  -  43) 
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Figure  3*3  HED  at  interface  3b. 


A^x33{p)  =  I  nfixp)  d$3  exp(-U3bz)  dk 
Ai33(p)  =  COS  <l>j  nf  \A,p)  exp(-U3bz)  dX 


(3-44) 


(3-45) 


The  boundary  conditions  to  solve  for  the  coefficients  of  (3  -  40)  -  (3  -  45)  are: 


@z  =  bn, 


(3-46) 


Plb  dz  P2b  dz 
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@  z  =  b2b 


Az13  _ 


lilb  |i2b 


1  C^A^13 

1  ' 

ElbPlb  9z 

e2bP2b 

9z 

Ax23  =  A^3 

1  ^Ax23 

1  ^Ax33 

_5(p) 

P2b  3z  lJ,3b  3z 

2jcp 

xb  xb 

Az23  _  ■^z33 

M-2b  M-3b 

1  3A^3 

1 

^A^3  _ 

e2bP2b  3z 

e3bP3b 

dz 

I  1 _ 1  \  ^^xl3 

UlbM-lb  ElbM-lbl 


’^x33 


(3  -  48) 


(3  -  49) 


(3  -  50) 


(3-51) 


(3-52) 


(3  -  53) 


and  the  Sommerfeld  radiation  condition  (3  -  29).  As  before,  the  term  5(p)/(2rtp)  in 
(3-51)  represents  the  unit  current  source,  this  time  at  interface  2b.  Solving  for  the 
unknown  coefficients  of  these  equations  yields 


H? t^p)  sinh(uibZ)  <1X 


(3  -  54) 


Azi3{p)  =  ^|  hH^p)  "f It cosh(uibz) dX 


4k 


D^eWo^X) 


(3-55) 


A"x23(P)  = 


4k 


H?\Xp) 


J  C 


Nbx3(x)  sinh(u2bz)  +  Ncx3(x)  sinh(u2bz) 

De(X) 


dX 


(3  -  56) 
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A^3(p)  = 


cos  <}> 

hP^Xd) 

Nbz3(^)  sinh(u2bz)  +  Nczsi^)  sinh(u2bz) 

471 

J 

c 

d^IxId^x) 

dk 


(3  -  57) 


Ax33(P)  = 


exp(-U3bz)  dX 


(3  -  58) 


b  COS(}> 

A^3(P)  =  -^ 

The  remaining  parameters  of  (3  -  54)  -  (3  -  59)  are  given  in  Appendix  C. 

3.2.3  HMD  at  interface  lb.  The  electromagnetic  fields  in  region  b  due  to  the 
equivalent  magnetic  current  on  the  ground  plane  (see  Figure  3  -  4)  are  found  with  (3  -  II) 
and  (3  -  12)  along  with  (3  -  6)  and  (3  -  8).  The  analysis  is  begun  by  finding  the  fields 


Hf\Xp) 


nU>-) 

d^xId^x) 


exp(-U3bz)  dX 


(3  -  59) 


z 


Figure  3  •  4  HMD  at  interface  lb  (dielectrics  1b  &  2b  are  transparent  for  clarity). 
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created  by  a  HMD  along  the  ;c-axis  at  the  ground  plane  with  a  unit  moment  Vdx  =  1  V  m. 

As  above,  the  x  and  z  components  of  the  electric  vector  potential  in  each  dielectric 
and  free  space  of  region  b  can  be  expressed 


Fxii{p)  =  sinh(uibz)  +  e$i  cosh(uibz)]  dX 

Fzi  i{p)  =  cos  (t)J  nf  \Xp)  a^i  sinh(uibz)  dX. 

Fx2i(p)  =  HH^p)[bxi  sinh(u2bz)  +  Cxi  cosh(u2bz)]  dX 

Pz2i(p)  ==  cos  <pj  Hf\xp)[bzi  sinh{u2bz)  +  Czi  cosh(u2bz)]  dX 


Px3l(p) 


=  |hJ?\xp) 


dxi  exp(*u3bz)  dX 


F^ i{p)  =  cos  <(>j  Hf  \xp)  dzi  exp(-U3bz)  dX 


(3  -  60) 

(3-61) 

(3-62) 

(3-63) 

(3-64) 

(3  -  65) 


The  boundary  conditions  to  solve  these  equations  for  the  unknown  coefficients 


are: 

N 

II 

o 

1  apj,,  _ 

5(P) 

eib  dz 

2jtp 

(3  -  66) 

@z  =  bn, 

^xll  -  *^x21 

(3  -  67) 
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1  9f5ii  . 

1  SfSzi 

Elb 

dz  e2b  dz 

(3- 

•68) 

=:£i2L 

Elb 

£2b 

(3- 

-69) 

1 

^F^^n 

.  1 

aF^2t  ■ 

^  e-^r\e*  A  /  1  1 

ElbM-lb  5z 

e2bP2b 

dz 

“  CUa  1 

\e2bP2b  EibM-lb 

1  ap 

(3- 

-70) 

@  z  =  b25 

Fx21=Fx31  (3-71) 

L  ^Fx21  _  1_^?3L 

eib  dz  esb  dz  (3  -  72) 

e2b  e3b  (3  -  73) 

_i_3eiL._L_?^  =  cos«|-J _ 1-1^ 

e2bli2b  e3bli3b  Ie3bli3b  e2b42b/  Op  (3-74) 

and  the  Sommerfeld  radiation  condition  (3  -  29)  with  A  replaced  with  F.  Similarly  as 
before,  the  term  -5(p)/(27tp)  in  (3  -  66)  represents  the  unit  magnetic  current  source  on  the 
ground  plane  (interface  lb).  The  solutions  of  (3  -  60)  -  (3  -  65)  are  then 


ri„(p)=^ 


sinKuikZ)  + 

u.i,  Dyx.) 


dX 


(3  -  75) 


Fzii{p) 


b^'Ub^'T  ^‘'’•^“ibz)  dX 

D^Xloax) 


(3-76) 
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Oft  p)  [Nb}ii(>-)  sinHuabZ)  ■■■  N?«iM  cosh(u2bz)] 

°  dSjIx) 


dX 


(3  -  77) 


f"2i(p)  = 


Hf\Xp) 


Nbzi  sinh(u2bz)  +  Nczi  sinh(u2bz) 


D^XlDyx) 


dX 


(3  -  78) 


Fx3i{p)  =  :tl 

d'J[x) 


4k 


exp(-u3bz)  dX 


(3  -  79) 


(3  -  80) 


The  remaining  numerator  terms  are  in  Appendix  C. 

32.4  HED  and  HMD  in  region  a.  The  solution  for  the  magnetic  vector  potential 
for  a  HED  at  interface  2a  is  much  simpler  since  there  is  only  one  dielectric  layer  in 
region  a.  The  vector  potential  for  the  single  layer  case  can  be  found  by  either  reducing 
(3  -  32)  -  (3  -  37)  or  (3  -  54)  -  (3  -  59)  by  letting  b2b  =  b^,  =  b^^,  U25  =  U3j,  =  02^.  e2b  = 


^3h  =  ^2a’  ^^b  =  ^^3b  =  ^^a>  “lb  =  “la’  ^Ib  =  ^la’  ^^Ib  =  ^^la’  o'" 

boundary  conditions  for  a  single  layer.  This  has  been  accomplished  by  other  researchers 

[5;  6;  12]  and  only  the  results  will  be  presented  here: 


Aii2(p)  = 


4jt 


H(?\^p)  sinh(uiaz)  dX 

dJx) 


(3-81) 


Ali2{p)  = 


d3(x)d*4x) 


COSh{uiaZ)dX 


(3  -  82) 
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Ax22{p)  =  ^1  Hq^^Xp)  exp(-U3bz)  dX 


A^2{p)  =  Hf ^Xp)  exp(-U3bz)  dX 

D^X)D“niX) 


(3  -  83: 


(3-84 


where 


Ula  =  V  X^  -  ,  U2a  =  V  X^  -  k^a 

with  kj^  and  k2^  defined  similarly  as  before  with  k2^  =  k^jj,  assumed  to  be  free  space 
k^.  The  solutions  for  the  HMD  at  interface  la  are 


Fxii(p)  =  ^  hHM  dX 

L  *  UiaDnJ(X/ 


(3-85 


f;u{p)  =  ^  HfW)-#¥^sinli(u,.z)ca 

DixpUxj 


F;2i(p)  =  3l|  H?>(Xp)?^|^ex[<.U3bZ)dX 


f;2.(p)  =  5^  HfW)^^exp(.„3bZ)dX 


(3-86 


(3-87 


with  the  zeros  of 


D^X)  =  Pal2U2a  +  UuCOth(biaUia) 


=  £al2  U2a  +  Uutanh(bi,Uu) 


(3-88 


(3-89 


(3-9C 
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defining  the  surface  wave  poles  in  region  a  where 


eai2=— ; 

£2a 


Ual2  =  — 
M-2a 


The  remaining  parameters  are  in  Appendix  C. 

3.2.5  The  scalar  potentials.  To  calculate  the  E  and  H  fields  of  (3  -  5)  and  (3  -  8), 
the  electric  and  magnetic  scalar  potentials  are  also  needed.  These  can  be  found  with  the 
Lorentz  gauge  conditions  for  electric  and  magnetic  sources: 

V- A  +  jo)ep.V  =  0  (3-91) 

V- F -»-jO3£^lVm  =  0  (3-92) 

When  applied  to  the  vector  potentials  calculated  above,  (3-91)  and  (3-92)  will  actually 
yield  the  potentials  of  two  opposite  point  charges;  electric  charges  for  (3-91)  and 
magnetic  charges  for  (3  -  92).  To  find  the  potential  Green’s  functions  of  (3  -  10)  and 
(3  -  12),  the  potentials  of  single  point  charges  are  needed.  The  potential  of  a  single 
electric  point  charge  and  of  a  single  magnetic  point  charge  are  found  using  the 
electrostatic  and  magnetostatic  relationship  linking  the  point  charge  and  dipole  potential: 


for  the  electric  point  charge  [12]  and 


(3  -  93) 


V  = 

(3  -  94) 

for  the  magnetic  point  charge.  By  integrating  V  or  over  x,  or  can  be  found 
within  an  arbitrary  constant.  Since  the  gradient  of  or  is  taken  to  calculate  the 
electric  or  magnetic  field,  the  value  of  the  constant  is  immaterial. 

3. 2. 5.1  Scalar  potentials  for  HED  at  interface  2b.  Applying  (3-91)  to 
(3  -  32)  and  (3  -  33)  yields 


which  results  in 


jweibHibV^2  =  - 


0Z 


(3-95) 
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=  HfW) 


- Uib - 


D^(x)  Dax)D^j[x). 


sinh(uibz)dX 


(3-96) 


As  stated  above,  this  is  the  potential  for  two  point  charges.  The  potential  for  a  single 
point  charge  is  found  by  applying  (3  -  93)  to  (3  -  96)  to  get 


v5i2(p)  = 


1 


4jtjcoeib|iib 


H?\Xp) 


uibN^iW 


Lo^eW  X  D^eWDy)^). 


sinh(uibz)dX 


(3  -  97) 


The  potentials  in  the  other  mediums  are  found  similarly: 


'V^22(P)  = 


± 


47CjO)e2bP2b 


H?’(Xp) 


n£«2(x) 

U2b  nUx) 


D^X)  X  D^xId^xU 
|,|nUx)  U2b  nUx) 


J  c 


1  L  D^e{x)  X  D^(x)Dyx)J 


sinh(u2bz) 

cosh(u2bz) 


dX 


^32(p)  = 


47tjO)e3bP3b 


H?\Xp) 


nSx2(x)  t  U3b  nSz2(x) 

.  D^x)  X  D^(x)Dyx). 


exp(-U3bz)  dX 


(3  -  98) 


(3  -  99) 


3. 2.5. 2  Scalar  potentials  for  HED  at  interface  3b. 


v5i3{p)  =  r^-‘^ - 1 

^  47tj{oeibHib  '  ^ 


N^3(x)  uib  N^X) 

Lo^elx)  X  Dax)D^Jx). 


47tj(0e2bP2b 


H^Xp) 


Nbx3(x)  U2b  Ncz3(x) 


D^x)  X  Dax)Dyx)J 

N^x3(x)  U2b  N^z3(x) 


sinh(uibz)dX 


sinh(u2bz) 


(3  -  100) 


D^(X)  X  D^xId^X). 


dX 


cosh{u2bz) 


(3-101) 
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v533(p)  = 

V^mqlllp) 

^^q2l{p) 

^^q3l{p) 

V"ql2(p)  = 

V‘q22(p)  = 

^^qll(P)  ■ 


4KjO)e3bP3b 


hH^p) 


NdxsIM  ^  Ulh  Ndxdi^) 

D^X)  X  caxlDUx). 


exp(-U3bz)  dX 


3. 2. 5. 3  Scalar  potentials  for  HMD  at  interface  lb 

^  feifr)  uid  nU>-) 

0  '  \  VTb  ^  \ 


(3  -  102) 


4KjcoeibPib 


cosh{uibz) 


f 


4:cjcoe2bP2b 


HoW) 


*  sinh(u[bz) 

dSIx) 

Nb,i(>.)  U2b  N^zi(x) 

.  D^x)  X  D^xId^X) 

^N^xi(x)  U2h  N^,i(x) 


dX 


/  L 

LdMx)  X  Dax)D‘J[x). 

I  nS,i(x)  ^  uxb  nSxA) 

[d^x)  X  Dax)Dyx). 


sinh{u2bz) 

COSh{u2bZ) 


- 1  hH>-p) 

t7Cj(l)e3bP3b  ' 

3.2 .5.4  Scalar  potentials  for  HED  at  interface  2a. 


exp(-U3bz)  dX 


(3  -  103) 


dX 


(3  - 104) 


(3  -  105) 


1 


47tjcoeiaP.i, 


X 


47tjO)e2#2a 


H?\Xp) 


H?^Xp) 


uu  NL2(x) 


L  d3(x)  X  d31x)d“4x)J 

Ndx2(M  ,  U2a  Ndz2(x) 

.  dJx)  X  d2(x)d“J(x). 


sinh(uiaz)dX 


exp(-U3bz)  dX 


(3  - 106) 


(3  - 107) 


32.53  Scalar  potentials  for  HMD  at  interface  la. 

f 

-1 - 1  H^W) 


4nj(DeiaHu 


N*xi(x)  uu  n;,i(x) 


D‘Jx)  X  d3(x)d*Jx)J 

+  sinh(uiaz) 

D^X) 


COSh(uiaZ) 


dX 


(3  - 108) 
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V?nq2l(p)  -  TTirdt 


47ljCO£2aU2a 


1  U2a  Ndzi(x) 


D^jx)  X  Dt(x)D?Jx)J 


exp(-U3bz)dX 


(3  -  109) 

All  vector  and  scalar  potentials  derived  above  were  checked  to  insure  they  met 
the  appropriate  boundary  conditions. 


Table  3  -  2  Necessary  Tangential  Fields 


HED  at  interface  2b  (Patch  1): 

HMD  at  interface  lb: 

1)  at  interface  1 

1)  at  interface  1 

2)  at  interface  2 

2)  E^^'^  at  interface  2 

3)  E^^  at  interface  3 

3)  E'^  at  interface  3 

1)  at  interface  1 

1)  at  interface  1 

2)  E^^”  at  interface  2 

2)  E^^  at  interface  a 

3)  E^^  at  interface  3 

HMD  at  interface  la: 

1)  at  interface  1 

2)  E^^"  at  interface  a. 

3.3  Constructing  the  Green’s  Functions 

From  the  vector  and  scalar  potentials  calculated  above,  the  E  and  H  fields  in  any 
medium  from  any  of  the  above  sources  are  found  using  (3  -  5)  to  (3  -  8).  However,  only 
the  tangential  fields  of  Table  3  -  2  are  needed  for  this  analysis.  In  the  following  sections, 
the  Green’s  functions  due  to  each  source  necessary  to  calculate  these  fields  are  construct- 
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3.3.1  E  fields  due  to  electric  sources  and  H  fields  due  to  magnetic  sources. 
Equations  (3  -  5)  and  (3  -  8)  are  used  to  calculate  the  E  and  H  fields  due  to  their 
respective  sources.  The  final  form  of  the  equations  is  independent  of  the  type  of  source 
(electric  or  magnetic)  because  of  duality  and  independent  of  the  source  and  observer 
interfaces.  The  Green’s  functions  for  the  case  of  an  electric  source  on  interface  2b  and 
observer  on  interface  3b  are  found  to  demonstrate  this  form. 

The  tangential  E  field  at  interface  3b  is  found  by  putting  (3  -  36)  and  (3  -  99)  into 
(3  -  5)  to  get 

E32"  =  Ay32y )  -  (z=bj^)  (3  - 1 10) 

»  b 

where  ^y32  is  the  magnetic  vector  potential  for  a  y  -directed  current  source  and  is  equal 
to  (3  -  36).  Since  only  a  unit  point  source  at  p'  =  0  is  assumed  (for  now),  it  is  seen  that 
the  integrals  of  (3  -  9)  and  (3  -  10)  over  the  point  source  leave 

G^2(P)  =  A$32(p)  (3-111) 

gS2(P)  =  a532(p)  (3-112) 

G$32(p)  =  V$32(p)  (3-113) 

where  G^2(P)  is  defined  as  the  magnetic  vector  potential  in  the  jc-direction  (first  x  in 
the  superscript)  on  interface  3b  for  an  infinitesimal  jr-directed  current  source  (second  x 
in  the  superscript)  on  interface  2b,  with  observer  at  radial  position  p.  The  Green’s 
functions  in  (3- 112)  and  (3- 113)  are  defined  similarly.  By  taking  the  gradient  in 
cylindrical  coordinates  and  using  the  identity  p  =  x  cos  (}>  +  y  sin  ()>,  (3  -  110)  becomes 

E5r(p)  =  cos  4,  T  +  -joG^^p)  -  sin  4,  y  ^ 

Notice  that  the  Green’s  functions  depend  only  on  the  radial  separation  p  =  |pl  between  the 
source  and  observer  and  not  the  relative  angular  positions.  This  expression  gives  the 
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tangential  E  field  anywhere  on  interface  3b  for  a  point  source  at  p  =  0  on  interface  2b. 
For  arbitrary  source  and  observer  points,  let  p  and  (])  become  (see  Figure  3-5) 


R=|p-pi 

„  .  ,  p  sin  <|)  -  p'  sin  d)' 

C  =  sin-*[^^ - ^ - 

Then  for  a  distribution  of  sources  over  interface  2b,  (3-114)  would  become 


(3-115) 


E32“(p)  =  f-j0)|  G^"2(R)J.2(p')ds'-cosc|  G5.32(R)q2(p')ds-  I 

+  Ga»2(R)Jy2(p')<is'-«nc|  0$.32(R)q2(p')ds'  y 


(3-116) 


where 


Qb  2  ._aGS32(R) 
Gq'32(R) 


By  properly  interchanging  the  Green’s  functions  and  source  terms  in  (3  - 116),  all  the 
tangential  E  fields  produced  by  HEDs  and  all  tangential  H  fields  produced  by  HMDs  for 
Table  3-2  can  be  found.  All  the  Green’s  functions  are  summarized  in  Appendix  D. 


Figure  3  -  5  Geometrical  relationships  for  x-directed  HED  at  arbitrary  p'  [i2]. 
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3.3.2  E  fields  due  to  magnetic  sources  and  H  fields  due  to  electric  sources.  The 
tangential  E  fields  at  the  dielectric  interfaces  due  to  the  magnetic  sources  on  the  ground 
plane  are  found  with  (3  -  6).  Again,  the  form  of  the  equations  are  the  same  for  either  E 
or  H  fields  so  only  an  example  is  done  here. 

The  tangential  E  fields  at  interface  2b  due  to  an  infinitesimal,  jc-directed  magnetic 
current  source  on  interface  lb  (ground  plane)  can  be  found  by  first  expressing  F21  in 
cylindrical  coordinates 


fS2i{P)  = 


cos 


4k 


-  (  sinh(u2bz)  +  NcxiM  cosh(u2bz)]  ^ 

Dy>.) 


(3-117) 


Fo^i(p)^  '^^”<l>[  sinh(u2bz)  +  Ncxi(A.)  cosh(u2bz)]  ^ 

^  D^X) 


=  HfW) 


(3-118) 


b  b 

Nbzi  sinh(u2bz)  +  Nczi  sinh(u2bz) 


D^(x)Dyx) 


dX 


(3-119) 


and  then  applying  (3  -  6)  in  cylindrical  coordinates  to  get 


El.T{p)  =  -7^ 


i 

e2b 


[IP  3$ 


b  \ 


3z 


-  .  (^Fp2i 

dz  dp  , 


P+! 


<> 


l(z=bib) 


(3  -  120) 


By  using  p  =  x  cos  <})  -b  y  sin  <j)  and  <|>  =  -x  sin  (J)  y  cos  <j>,  the  x  and  y  components  of  the 
E  field  can  be  found: 
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E5r(p)= 


sin  20 


X.  Nbzi  sinh(u2bbib)  +  Nczi  sinh(u2bbib) 

47tE2b^  /  b  b  I 

.1  H^2)j^p)Nbzi  sinh(u2bbib)  +  Nczi  sinh(u2bbib) 

P  ^  d^IxId^x) 


(3-121) 


H^^\Xp)  u^^^^bxi(x)  cosh(u2bbib)  +  Ncxi(x)  sinh(u2bbib) 

D^X) 

^<bl  sinh(u2bbib)  +  Nczi  sinh(u2bbib)  \ 

Daxwx) 


-  cos 


^cos^  jNbii  smli(u2bbib)  ■^  N^i  sinKuibbib) 

j  ‘  d^xIdUx) 

J  C 


(3  -  122) 


where  E21*  denotes  a  ^-directed  field  for  an  x-directed  magnetic  current  element.  Notice 
again,  the  integrals  of  (3  -  121)  and  (3  - 122)  depend  only  on  the  radial  distance  between 
the  source  and  observer  points.  For  an  arbitrary  distribution  of  x-directed  magnetic 
current  elements  and  using  (3  - 115),  (3  - 121)  and  (3  - 122)  are  represented  as 


|Gg;(R’C)'<i(p')<*s' 


(3  -  123) 


(3  -  124) 


where  -Ge2i(R.C)  equals  the  right-hand-side  of  (3-121)  and  -GgjjlR.C)  equals  the 
right-hand-side  of  (3  - 122)  for  arbitrary  source  and  observer  positions. 

The  E  fields  for  y-directed  magnetic  current  sources  are  found  by  letting  0 
become  0  -  n/2  in  (3  -  1 17)  to  (3  -  1 19)  [I2]  and  repeating  the  same  procedure  as  above. 
The  X  and  y  components  of  the  E  field  at  interface  2b  due  to  an  infinitesimal,  y  -directed 
magnetic  current  source  on  interface  lb  are  then 
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pbyy  _ 
^21  “ 


.i(  n)?.  sinh(u2bbib) -^  N^,1  sinh(u2bbib)  ^ 

sin  2$)  "  dKxWJi) 

4K£2b  f 

+  J-  H^%p)^bzi  sinh(u2bbib)  +  N^zi  sinh(u2bbib)  ^ 

P  ^  D|(XW\) 


cosh(u2bbib)  +  N^,i(x)  sinh(u2bbib)  ^ ' 


(3  -  125) 


fcos^J  Hi,MXp)x  sinh(u2bb^b)-b  sinh(u2bbib)  ^ 

)  cos  2<{>  ^(2)|n  sinh(u2bbib)  +  N^^i  sinh(u2bbib)  ^ 

P  ^  D^>.WX) 


(3  - 126) 


which  are  the  same  as  (3  - 121)  and  (3  - 122)  except  for  a  few  sign  changes.  For  an 
arbitrary  distribution  of  y  -directed  magnetic  current  elements,  the  E  fields  are  represent- 


Ejr(p)=- [Ggi[(R,C)M5,(p')ds' 

^;'’(P)=- foSjlR’ClMSlIpIds' 


(3  -  127) 


(3  - 128) 


The  total  tangential  field  can  be  written  more  compactly  as 


where 


E2'"(p)=-|G|2,(R,C).M5(p')<is' 

5|2.(R.C)4l?-^l!! 

LGE2i(R.C)yx  GE2i(R,C)yy 


(3  -  129) 


(3  - 130) 
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The  tangential  E  fields  at  interface  3b  for  the  magnetic  currents  are  found  by 
replacing  bj^j  with  b2|j  (3-1 17)  to  (3  -  1 19)  and  repeating  the  above  procedure  to  get 


and  in  region  a 


EsHp)  =  - 1  fi3i(R.C)-  M^(p')  ds' 
E2r{p)  =  -|^2i(R.C)-Ml(p')ds' 


-  131) 


(3  -  132) 

Starting  with  (3  -  7)  and  following  the  above  procedure,  the  tangential  H  fields  at 
the  ground  plane  are  found  for  electric  currents  on  the  dielectric  interfaces.  For  currents 
on  interface  2b 


currents  on  interface  3b 


H?f{p)  =  |Sii2{R.C}  J2(p')ds' 
H?np)  =  |Sii3(R.i;)-  J3(p')ds' 


(3-133) 


(3  -  134) 


and  the  scattered  currents  on  interface  2a 


H!r(p)=  ^i2(R,C)-  Jf(p')ds' 

h  (3  - 135) 

This  last  equation  can  also  be  used  to  find  Hi2"’(p)  for  the  incident  current,  The 

complete  equations  for  the  Green’s  functions  of  (3  - 129)  to  (3  -  135)  are  in  Appendix  D. 


3.4  Matrix  Equations 

The  unknown  current  and  charge  distributions  on  the  conductors  and  aperture  are 
found  by  enforcing  the  boundary  conditions  of  Table  3-1.  The  boundary  conditions  are 
expressed  as  integral  equations  involving  the  calculated  Green’s  functions  and  unknown 
distributions.  Using  appropriate  basis  and  testing  functions,  the  integral  equations  are 
converted  into  matrix  equations  that  can  be  solved  for  the  unknowns. 
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3.4.1  The  integral  equations.  To  meet  the  boundary  condition  =  0  on  patch 
1,  the  sum  of  the  tangential  electric  fields  due  to  the  sources  on  patch  1,  patch  2  and  the 
equivalent  sources  on  the  aperture  must  equal  zero.  Using  the  Green’s  functions  defined 
earlier  and  (3  -  5),  this  is  expressed  as 


jco  ^2(plp')- J2(p')cls'-Vj  G522(plp')q2(p1ds' 

Jpl  Jpl 

)|  ^23(plp')- J3(p')ds'-v1  GS23(plp')q3(p')ds' 

Jp2  Jpl 

+  ^2i(plp')-Mi(p')ds'  =  0 

Jad 


-JCO 

-jco 


where 


GA22(plp')  = 


G^A^P')xx  0 

0  GA^22(p>p')yy. 

GA23(plpi  is  defined  similarly,  and  from  Figure  3  -  1 

mV)  =  -Mi(p') 

Similarly  on  patch  2 


-jcol  GA32(plp')-  J2(p')ds'- V  I  Gq32(plp')  qzip')  ds' 
Jn  Jpi 

-jcol  GA33(plp')- Jslp'lds'- v‘1  G$33(plp')  q3(p')  ds' 
JP2  Jf2 

+  1  ^3i(plp')-Mi(p')ds'  =  0 

Jap 


(3  - 136) 


(3  - 137) 


(3  - 138) 

On  the  feedline,  the  difference  between  the  tangential  E  field  due  to  the  incident 
current  and  the  sum  of  the  fields  due  to  the  reflected  current  and  equivalent  magnetic 
current  must  equal  zero: 
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if>|  ^2(plp')-  Jinc(p')  ds'-  V‘  I  Gq22(plp1  qinc(p')  ds' 

/feed  /feed 

+jco|  GA22(plp')-  Jf(p')ds'+ V  I  Gq22(plp')qf(p')ds' 

^foed  /iced 

+  |  G^2i(plp')-M,(p')ds'  =  0 
Ap 


(3-139) 


where 

M?(p')  =  Mi(p') 

Truncating  the  microstrip  feed  at  2  to  3  wavelengths  from  the  aperture  is  sufficient  [5]  to 
produce  accurate  results  for  the  type  of  basis  and  test  functions  discussed  below. 

The  final  boundary  condition  that  must  be  enforced  is  continuous  through 
the  aperture.  The  tangential  magnetic  field  at  the  aperture  due  to  all  sources  in  region  a 
must  equal  the  tangential  field  due  to  all  sources  in  region  b.  Using  (3  -  8)  and  the 
defined  Green’s  functions,  the  integral  equation  is 


I  GHnlpIp')- Jinc(p')ds'+ I  GHi2(plp')- Jf(p')  ds' 

/bed  /feed 

■jcoj  0?ii(p(p')  Mi(p')ds'-v‘j  G^„(plp')q„,i(p')ds' 

Jkp 

=  I  Gnulplp')-  Ji(p')ds'+  I  GHoipIp')-  J2(p')ds' 

/pi  /p2 


-jtoj  Sii(plp')-Mi{p')ds'+V‘j  G^n(plp')q^i(p')ds' 

jAp 


(3  -  140) 

3.4.2  Basis  and  test  functions.  To  solve  for  the  unknown  currents  and  charge 
distributions  in  (3  - 136)  to  (3  -  140),  a  method  of  moments  approach  using  subsectional 
basis  functions  [11]  can  be  applied  to  convert  the  integral  equations  to  matrix  equations. 
The  source  surface  (antenna  patch,  feedline,  or  aperture)  is  divided  into  elementary 
domains  (cells  Sj)  and  a  basis  function  is  defined  over  each  cell  [5].  The  observer  surface 
is  also  divided  into  cells  with  a  test  function  applied  to  each  observer  cell  Cj.  The 
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number  of  cells  the  source  or  observer  planes  are  divided  into  depends  on  the  area:  a 
small  surface  will  take  fewer  cells  to  describe  than  a  larger  surface.  Hence,  the  size  and 
number  of  cells  used  to  divide  the  source  and  observer  planes  may  be  different. 

Mosig  showed  in  [5]  that  overlapping  rooftop  functions  can  be  used  effectively  to 
expand  the  x  and  y  components  of  the  surface  currents.  The  basis  functions  for  the 
surface  charges  are  then  2-D  pulse  doublets,  according  to  the  continuity  equations 
(see  Figure  3-6).  A  1-D  pulse  test  function  can  be  used  to  greatly  simplify  the 
calculation  of  the  matrix  elements  without  sacrificing  convergence  time  or  accuracy. 
This  choice  also  eliminates  the  need  to  compute  field  values  near  the  surface  edges, 
where  field  singularities  can  adversely  affect  the  performance  of  the  moment  method  [13]. 


1  Basis  Functions  | 

Test 

Function 

Current 

Charge 

•a  N 

-a  1  +a 

2D-Rooftop 

-a/2  1  +a/2 

1-D  Pulse 

Figure  3  -  6  Basis  and  test  functions.  The  two-dimensional  functions  are  independent  of  the 

transverse  coordinate  [sj. 


The  vector  rooftop  function  associated  with  each  current  source  cell  Sj  is  defined 

as  Tj.  To  adequately  describe  the  current  on  a  conductor  or  over  the  aperture  requires  N 

x-directed  functions  and  N„  y-directed  functions  for  a  total  N  =  N  +  N  .  Therefore, 

y  ^  y 


-42- 


III.  Theory 


Tip')=|y 
\0; 


1-! 


!-'• 


(p'  -  P'j)l 
'■  (p'  -  P'j)l 


;  |x- (p'-p'j)|<a.|y-(p'-p'j)|<b/2,  l<j<Nx 

;  k-(p'-p'j)|<a/2,|y-(p'-p'j)|<b,  Nx<j<N 


(3-141) 


Otherwise 

The  X  and  y  components  of  the  current  are  expanded  as 

Js.  =  ffiajTj(p')  Js,  =  i  i  «j Tj(p') 

°j=i  j=N.+i  (3  -  142) 

where  ttj  are  the  unknown  coefficients,  p  j  is  the  center  of  the  j  current  cell  on  the 
source  plane,  b  is  the  dimension  of  the  charge  cell  perpendicular  to  the  x  -direction,  and 
a  is  the  dimension  perpendicular  to  the  y  -direction  (see  Figure  3  -  7a)  [13].  Note  in 
Figure  3  -  7a,  the  jc-current  cell  Sj  and  y  -current  share  the  same  charge  cell.  In 

fact,  any  charge  cell  can  belong  to  as  many  as  4  current  cells:  2  x-current  and  2  y  -current 
cells.  Thus,  the  current  cells  automatically  overlap.  The  continuity  equations  can  be 
used  to  expand  the  surface  charge  density,  yielding 


N 


“irtp'l 


(3  -  143) 


where  n/p') = -  V  which  is  a  two  dimensional  pulse  doublet  function  (see 
Figure  3  -  7b)  defined  as 


+  1;  0  <  X-  (p'  -  p'j)  <  a,  I y •  (p'  -  p'j)l  <  b/2,  1  <  j  < 

-  1;  -  a  <  X-  (p'  -  p'j) <  0,  |y-  (p'  -  p'j)l  <  b/2,  1  ^  j  ^  N* 

nj(p')  =  \  +  1;  lx-  (p'  -  p'j)l  <  a/2, 0  <  y-  (p'  -  p'j)  <  b,  N*  <  j  <  N 

-  1;  |x-(p'-p'j)|<a/2,-b<y-(p'-p'j)<0,  Nx<j^N 

\0;  Otheiwise  (3  -  144) 

These  functions  are  used  similarly  to  expand  the  equivalent  magnetic  currents  and 
charges  of  the  aperture. 

The  relationship  between  the  number  of  charge  and  current  cells  is  not  exactly 
simple,  since  the  relationship  depends  on  the  shape  of  the  conductor  or  aperture.  For  a 
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x-current  cell 


Figure  3  -  7a  Segmentation  of  source  plane  into  elementary  charge  and  current  cells. 


Figure  3  •  7b  x-directed  current  cell  at  p^O  and  associated  surface  charge  [13]. 

rectangular  patch  or  aperture  with  m  xn  charge  cells,  the  number  of  x  -directed  current 

cells  is  N„  =  (m  -  l)n  and  the  number  of  y-diiected  current  cells  is  =  m(n  - 1)  [13]. 

^  y 

3.4.3  Convert  integral  equations  to  matrix  equations.  The  basis  and  test 
functions  defined  above  are  applied  to  the  integral  equations  (3  - 136)  to  (3  - 140)  to 
define  a  series  of  linear  equations  that  can  be  solved  for  the  unknown  coefficients  using 
standard  matrix  techniques.  The  test  function  is  applied  to  the  current  cel..'?  of  the 
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observer  plane  (unprimed  coordinates)  and  the  basis  functions  are  applied  to  the  current 
cells  of  the  source  plane  (primed  coordinates).  The  integral  equation  (3  - 136)  then 


becomes 


il[r[  d*-  [  ^2(plp')-T2j(p')ds' 

i=i  j=i  k 

jcoa^^l  [^S22(p2ilp')  -  G522(p2ilp')]  n2j(p')  ds' 

Js2i 

iii'rf  <“■  f  0«3(plp'lT3ip-)ds' 

i-1  i-1  [  ’Vc,  A, 

+  I  [^S23(p2ilp')  -  GS23(p2ilp')]  risjlp')  ds' 

Js3i 

-  i  i  [rf  dl-  (  S2i(plp')-  Tij(p')ds'ja|  =  0 

i=l  j=l  A.i 


(3  - 145) 


The  index  variable  j  always  refers  to  elements  on  the  source  plane  and  the  index  variable 
i  always  refers  to  elements  on  the  observer  plane.  The  number  N2  is  the  total  number  of 
current  elements  (,x  and  y-directed)  on  the  conductor  at  interface  2b  (patch  1)  and  Nj  and 
arc  defined  similarly.  The  other  variables  are  defined 


=  forl^j<N2x 
\a2  forN2x+l^j^N2 

C2j  =  P2j  to  P2j 

with  tjj,  t^j  and  Cij,  C^j  defined  similarly. 

The  matrix  form  of  (3  -  145)  is 

[4“lk^]-[z5’’][ar]-[cff][aj]  =  0  (3-. 46) 

where  the  [Z]  components  are  impedance  matrices  in  units  of  ohms/length  and  the  [C] 
component  is  a  coupling  matrix  in  units  of  1/lcngth.  The  column  vectors  and 
are  the  amplitudes  in  units  of  amps,  of  the  electric  current  elements  on  the  conductors  at 
interface  2b  (patch  1)  and  3b  (patch  2),  respectively.  The  column  vector  [«j]  represents 
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the  amplitudes  of  the  equivalent  magnetic  current  elements  of  the  aperture  in  units  of 


volts.  The  X  and  y  components  of  the  matrices  and  currents  are  related  by 


^xx  i^xy 
^ryx  ^yy 


a, 


1“^ 


=  or 


where  K  is  any  impedance,  admittance,  or  coupling  matrix.  The  individual  components 
of  the  matrix  are  determined  by 


Kf 


v-y* 


for 


for 


'  (N,<j<N 

'  \N*<j<N 


(Nx<i^N 

ll^j<Nx 

a*  for  1  <  j  <  Nx 

aj  for  Nx  <  j  ^  N 

In  this  way,  the  x  and  y  components  of  the  field  for  all  x-  and  y  -directed  components  of 
the  surface  current  are  described  with  one  matrix  operation. 

The  matrix  forms  for  (3  -  138)  to  (3  - 140)  arc  then 

[zr/l[ari-[4i”ik’]-[cr/'][aj]=o  o.mt) 

[C^'l[a]].[z^fa3  =  [z^^][arl  (3-148) 

[c5'fan>k‘far]^|[<']  -[Vtj"]!k]-[cf Jk]  =  [c,f][«rl  (3-149) 


where  the  [Y]  components  are  admittance  matrices  in  terms  of  mhos/length.  The  column 
vector  [®j]  are  the  amplitudes  of  the  reflected  electric  current  elements  on  the  feedline 
and  [®j"i  are  the  amplitudes  of  the  known,  incident  electric  current  elements.  The 
equations  for  the  various  impedance,  admittance,  and  coupling  matrix  elements  are 
summarized  in  Appendix  E. 

3.4.4  Solution  matrix.  The  matrices  defined  in  (3  - 146)  to  (3  - 149)  can  now  be 
used  as  sub-matrices  in  the  solution  matrix 
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vl/1 
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“j. 

[cfj'"] 

[0] 

[0] 

( 

cf  ]  [z,f2]  . 

[af] 

(3  -  150) 


where  [yJ,- J  “  [y?/  +  [Y^/  j.  All  the  diagonal  sub-matrices  represent  tangential  fields 
created  by  sources  in  the  same  plane.  Each  diagonal  term  of  these  sub-matrices 
determines  the  tangential  field  at  that  segment  created  by  the  source  in  that  segment. 
Therefore,  the  diagonal  terms  of  the  entire  matrix  represent  the  tangential  fields  at  every 
segment  due  to  the  sources  in  that  segment.  Since  these  fields  are  the  strongest 
calculated,  the  matrix  is  diagonally  dominant.  The  off-diagonal  sub-matrices  represent 
the  tangential  fields  at  one  plane  created  by  sources  on  a  different  plane. 

The  observed  fields  for  the  various  source  planes  are  determined  by  the  row  and 
column  positions  of  the  matrix.  The  row  positions  determine  the  the  observed  field  and 
the  column  positions  determine  the  sources  of  the  field.  For  instance,  the  first  row  of 
sub-matrices  determines  the  tangential  electric  field  at  interface  2b  over  patch  1. 
Columns  1,  2,  and  3  designate  the  sources  of  the  field  as  those  on  interface  2b,  interface 
3b,  and  interface  1  (ground  plane),  respectively.  Note  that  sources  on  the  feedline  make 
no  contribution  to  the  field  at  interface  2b,  so  the  fourth  column  has  a  zero  sub-matrix. 

The  number  of  rows  and  columns  of  each  sub-matrix  are  determined  by  the 
number  of  segments  (current  cells)  the  observer  plane  and  source  plane  are  divided  into. 
The  number  of  columns  (/'  index)  will  equal  the  total  number  of  current  cells  on  the 
source  plane.  The  number  of  rows  (i  index)  equals  the  total  number  of  current  cells  on 
the  observer  plane.  Therefore,  the  sub-matrices  in  each  row  of  the  solution  matrix  have 
the  same  number  of  rows,  but  a  different  number  of  columns.  The  diagonal  sub- 
matricies  are  square,  since  the  number  of  observer  and  source  cells  are  the  same.  This 
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results  in  a  square  solution  matrix  with  the  number  of  rows  or  columns  equal  to  the  total 
number  of  rows  or  columns  of  the  diagonal  sub-matrices. 


A^.  Application  and  Results 


To  apply  the  moments  method  procedure  outlined  in  chapter  in  and  determine  the 
surface  current  and  charge  distributions  of  the  antenna,  numerical  evaluations  of  the 
Green’s  functions  must  be  developed.  In  this  chapter,  the  functions  to  be  integrated  are 
first  investigated  and  characterized  so  proper  numerical  techniques  can  be  applied.  The 
actual  numerical  computations  are  done  with  custom  Fortran  code  combined  with  several 
commercial  subroutines  from  the  IMSL  libraries  [14].  These  commercial  routines  were 
used  because  of  their  proven  accuracy  and  to  save  time  in  code  development. 

4.1  Characteristics  of  the  Integrands. 

Because  of  the  complexities  of  the  integrands,  the  numerical  integrations  in  the 
Green’s  functions  are  far  from  simple.  The  integrands  are  defined  over  an  infinite  path  in 
the  complex  plane.  The  integrands  are  oscillatory,  they  decay  at  different  rates,  and  they 
involve  branch  cuts  and  possible  singularities  on  the  real  axis.  By  examining  the  various 
characteristics  of  the  integrands,  simplifications  are  made  to  aid  the  numerical  integration 
process.  Also,  the  potential  difficult  points  for  numerical  integration  are  located  and 
dealt  with. 

4.1.1  Deformation  of  the  integration  interval  into  the  real  axis.  The  integrals  of 
the  Green’s  functions  in  Appendix  D  can  be  written 


Since  /  is  always  an  even  function,  the  integration  path  C  can  be  deformed  into  the  real 
axis  and  (4-1)  becomes 
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Sn{f)  =  2j 

Jo  (4  -  2) 

The  function  S^(f)l2  is  then  the  Hankel  transform  of  X”/  [12].  If  the  function  /  has  a 
singularity,  the  integration  path  cannot  be  completely  deformed  into  the  real  axis.  The 
path  must  go  around  the  pole  on  the  half-plane  Im(X)  >  0.  Equation  (4  -  2)  is  then 
changed  to 

Sn(f)  =  2Pv[  J„(Xp)x"^M?^1ciX-2dS  RiJn(>.pip)C' 

Jo  ‘  (4  -  3) 

where  PV  denotes  the  principle  Cauchy  value  of  the  integral  [12]  and  Xp  -  is  the  location 

th 

of  the  pole.  The  residue  is  given  by 

Ri  =  Lim  (X  -  Xpi)  f(x  ) 

(4 . 4) 

4.1.2  Location  of  the  poles.  The  location  of  the  poles  in  the  integrands  are  deter¬ 
mined  by  the  zeros  of  the  denominator  functions  Dm  and  De  in  region  b,  and  Dm  and  D® 
in  region  a.  The  zeros  of  these  functions  are  located  on  the  real  axis  for  real  (lossless)  di¬ 
electric  permittivities  and  permeabilities.  Standard  numerical  techniques  are  then  used  to 
locate  the  zeros. 

The  location  of  the  zeros  are  further  restricted  through  closer  examination  of  the 
denominator  functions.  Looking  at  Dm  along  the  real  axis,  the  following  expressions  are 
obtained: 

for0<X<k3j,<k2i,<kiij; 

D^x)  =  ju3bl  -luibl  tar(bib  juibl  )]iu2bl  cos(lu2bl  tb) 

+  [-ebi2 1 U2bP  -  eb23 1 Uibl  j U3bl  tan(bib  I  uibl )]  j  sin(  |  U2bl  tb)  (4  -  5) 

for0<k3^<X<k2b<kib; 
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+  [-ebi2 1  U2bl^  -  Eb23 1  Uib!  U3b  tan(bib  I  uibl )]  j  sin{  j  U2bl  tb) 


(4-6) 


forO<k3j,<k2b<X<ki^,; 

D^j^^j^t^bB  U3b  -luibl  tan(bib|uib! )]  U2b  cosh(u2b  tb) 

+  [ebi2  uib  -  eb23  luibl  U3b  tan(bib  |  uibl )]  sinh(u2b  tb) 


(4-7) 


and  for  0  <  k^^^  <  k2jj  <  k^j^  <  X; 


D^X)  =  t^‘’i3U3b  +  uib  tanh(bibUib)]  U2b  cosh(u2b  tb) 

+  kbnuib  +  eb23UibU3b  tanh(bibuib)]  sinh(u2b  tb) 

If  k2jj  and  k^^j  are  reversed  (k^j^  <  k2jj),  only  (4  -  7)  needs  to  be  changed 


(4-8) 


j^yj^)_[ebi3U3b  +  uib  tanh(bibUib)]ju2blcos(ju2b|tb) 

+  kbBuib  +  eb23UibU3b  tanh(bibUib)]  j  sin(  | 


U2b|tb) 


(4-9) 


where  tjj  =  b25  *  bjb-  No  zeros  can  exist  in  (4  -  5)  since  both  the  real  and  imaginary 
parts  of  the  equation  would  have  to  go  to  zero  for  the  same  real  value  of  X.  There  are  no 


zeros  in  (4  -  8)  either,  since  the  function  is  always  positive  over  the  specified  interval  for 
X.  Zeros  can  exist  in  (4  -  6),  (4  -  7),  or  (4  -  9)  because  these  equations  are  all  real  or 
imaginary  over  the  specified  interval  and  contain  oscillatory  functions.  Therefore,  zeros 
can  only  exist  in  the  interval  k^jj  <  X  <  k^^^  where  k^^^^  is  the  maximum  of  k2jj  or  kjj^ 
(see  Figure  4-1).  Notice  that  the  real  part  of  Dm  goes  to  zero  at  X  =  k^j^,  and  while  the 
imaginary  part  is  continuous,  the  first  derivative  is  not.  Special  care  must  be  taken  when 
integrating  near  this  point.  It  appears  in  Figure  4  -  1  that  X  =  k2jj  is  a  possible  zero  of 
However,  U2f,  is  also  present  in  the  numerator  of  all  the  Green’s  functions  in  such  a  way 
that  the  integrand  is  always  totally  smooth  at  X  =  k25- 

It  has  been  shown  by  other  researchers  [9]  that  there  is  always  at  least  one  zero  for 
Dm.  In  other  words,  there  is  no  cut-off  frequency  for  the  TM^  surface  wave  mode.  For  a 
given  frequency  and  material  parameters,  Dm  =  0  can  be  solved  numerically  for  the  pole 
value  X^u.  The  FindRoot  function  in  the  commercial  software  package  Mathematica 
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Figure  4  - 1  Normalized  values  of  D^ksb  on  the  real  axis  for  (a)  f »  4  GHz,  =  1 .6  mm, 
b2b  =  ^  S  'Tim,  Jj  .  SEq,  Egb  =  2.5Eq,  Egj^  =  Eq,  4i  b  “  H2b  “  ^^b  =  ^'o■  (3) 

except  E^j,.2.5EQ,£2b-5Eo. 


[15]  and  the  Fortran  subroutine  dzbren  from  the  IMSL  libraries,  provided  very  accurate 
results  for  the  poles  of  For  example,  using  Mathematica  the  normalized  root  for 
Figure  4  -  la  is 

Xpj/k3j,  =  1.040360728 
Putting  this  back  into  yields 

D^3b  =  j  6.995646  X  10 
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The  normalized  root  for  Figure  4  -  lb  is 
\b^3h  =  1-057332766 

which  yields 

D^3b  =j  5.391632  X  10'^^ 

when  put  back  into  D^.  The  accurate  location  of  the  pole  is  critical  for  the  successful  nu¬ 
merical  integration  of  the  functions. 

There  may  or  may  not  be  a  zero  for  De  depending  on  the  observation  frequency. 
An  expansion  of  De  similar  to  (4  -  5)  -  (4  -  9)  shows  that  like  D^,  any  possible  zeros  must 
exist  on  k^jj  <X<  Figure  4-2  shows  a  normalized  plot  of  Dg  for  the  same  param¬ 
eters  as  Figure  4-1.  At  this  observation  frequency,  there  are  no  zeros  for  De,  except  at  X, 
=  k2jj;  but  this  does  not  produce  a  singularity  in  the  integrand  for  the  same  reasons  stated 
above  for  Dm-  It  was  shown  in  [9]  that  when  a  given  surface  wave  turns  on  =  k^jj. 
Making  this  substitution  into  Dg  yields 

(cDcesbjxsb  -  o^EibHiby^^  (cOcEsbUsb  -  0)ce2bli2by^^ 

*  COth(bib(o)ce3bli3b  -  Wceiblllb/^^)cosh(tb(c0ce3bli3b  -  t0cE2bll2bf^^) 

+  (wcEsbM-Sb  -  0)ce2bli2b)  Sinh(tb(c0ce3blt3b  -  0)ce2bM-2b)'^^)  =  0  (4  - 10) 

where  O)^  is  the  cut-off  frequency  of  the  TE  surface  wave  modes  for  the  given  material 
parameters.  Equation  (4  -  10)  is  a  transcendental  equation  that  can  be  solved  with  nu¬ 
merical  techniques  similar  to  those  used  to  solve  the  roots  of  D^.  For  example,  the  first 
cut-off  frequency  for  the  case  of  Figure  4  -  2a  is  /^  =  12. 1 138  GHz  and  for  Figure  4  -  2b 
is  =  7.93764  GHz. 

The  surface  wave  modes  are  similar  in  behavior  to  those  of  a  single  dielectric 
layer  [9].  The  dominant  TM  ^  mode  has  no  cut-off  frequency  and  the  following  modes  al¬ 
ternate  TEj,  TM2,  TE2,  etc.,  with  each  mode  turning  on  at  a  higher  frequency.  If  the  de¬ 
sign  frequency  of  the  antenna  is  kept  below  the  cut-off  frequency  of  the  first  TE  surface 
wave  mode,  only  the  one  singularity  of  D^  has  to  be  dealt  with  for  the  Green’s  functions 
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Figure  4-2  Normalized  values  of  De/ksb  on  the  real  axis  lor  (a)  f  =  4  GHz,  =  1 .6  mm, 
bgb  =  4.8  mm,  =  Se^,  -  2.56^.  £3^,  -  e^.  =  Hgb  =  ^'sb  *  ^^o= 

except  e^jj  =  2.5eQ,  e2b’=^^o 

in  region  b.  This  case  is  assumed  for  the  rest  of  the  analysis. 

The  analysis  of  the  poles  for  and  D|  is  covered  in  detail  in  [12],  yielding 
similar  results  as  above.  One  difference  is  that  for  the  single  layer  case,  the  cut-off  fre¬ 
quency  for  the  first  TE  surface  wave  mode  can  be  found  analytically  by 

f^  = - 1 - 

4bu((eai2-I)Eolio)’^^  (4-11) 
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For  the  same  first  dielectric  layer  used  in  region  b  in  Figure  4  -  la  (bj^  =  1.6  mm,  £319  = 
5)  the  cut-off  frequency  is  23.3955  GHz.  It  is  seen  that  a  proper  choice  of  material 
parameters  and  operating  frequency  avoids  operation  above  the  first  TE  surface  wave 
mode  in  either  region  and  greatly  simplifies  the  calculations. 

4. 1.2.1  Real  and  imaginary  parts  of  the  integrands.  Combining  the  anal¬ 
ysis  of  the  denominator  functions  as  shown  in  (4  -  5)  through  (4  -  9)  with  a  similar  analy¬ 
sis  of  the  numerators  of  the  integrands,  reveals  the  integrands  are  always  real  for  X  > 
in  region  b  and  X  >  ^2^  in  region  a.  The  integrands  axe  always  complex  (both  real  and 
imaginary  parts  present)  for  X  <  and  X  <  k2^-  Numerical  integration  routines  can 
thus  be  focused  efficiently  on  either  the  real  or  imaginary  components  present  in  the  dif¬ 
ferent  integration  intervals. 

4.1 .3  Rate  of  decay  of  the  integrands.  Since  all  the  integrands  for  the  Green’s 
functions  in  Appendix  D  oscillate  and  decay  as  X  00  (for  R  >  0),  it  can  be  shown  that 
the  integrals  will  converge.  As  X  ->  <»,  the  integrands  oscillate  about  the  X  axis  due  to 
the  Bessel  functions  and  decay  at  a  rate  that  is  determined  later.  If  a  function  is  integrat¬ 
ed  over  each  half-period  separately  as  shown  in  Figure  4  -  3,  and  these  values  summed  as 
a  series,  the  integral  is  represented  as 


S  =  I  (-ir'an 

n=l 


(4-12) 


where  a^  are  the  absolute  values  of  the  integrals  over  each  half-period.  Since  — >  0  as 
n  the  Leibnitz  criterion  (an  alternating  series  converges  only  if  the  terms  are 

monotonic  decreasing)  is  met  and  the  series  converges  [16]. 

The  rate  at  which  the  series  converges  is  dependent  on  the  rate  of  decay  of  the  in¬ 
tegrand.  The  integrands  for  the  different  Green’s  functions  decay  in  two  different  ways. 
Integrands  for  Green’s  functions  where  the  source  and  observer  are  in  the  same  plane 
decay  as  as  X  — >  <».  Take  for  example 
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Figure  4-3  Integration  over  discrete  haK-periods  of  oscillating,  decaying  function. 


r 


Jo(Xr)  X  [U2b  cosh(u2b  (b2b  -  bib))  +  ^^b23  U3b  sinh(u2b  (b2b  •  bib))]  ^ 
([^bi3U3b  +  Uib  coth(bibUib)]  U2b  cosh(u2b(b2b  -  bib))  I 
I  +  [l^bnuib  +  |ib23UibU3b  coth(bibUib)]  sinh(u2b(b2b  -  bib))l 


(4-13) 

which  represents  the  magnetic  vector  potential  on  interface  2b  at  a  distance  R  from  a 
point  source  on  interface  2b.  As  X,  — >  <»  the  integrand  becomes 


/nr  cos(XR  -  S-) _ [X  exp(Xtb)  -b  \ih23  ^  exF(>.tb)] _ 

{[nbi3  ^  +  A,  ]  X,  exp(Xtb)  +  .P-bnX  +  )ib23X  .  ex  p(Xtb))  (4  -  14) 


where  the  Bessel  function  has  been  replaced  by  its  asymptotic  approximation  [17].  This 
further  reduces  to 

^  cos(XR  -  S’) - ]  ^  - 

V  JtXR  '  ^^bl3  +  1  +  ltbl2  +  ltb23  (4  -  15) 

and  it  is  apparent  the  integrand  decays  algebraically  as  for  large  X.  This  is  true  for 

all  Green’s  functions  where  the  source  and  observer  are  on  the  same  plane. 
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The  integrands  for  the  Green’s  functions  where  the  source  and  observer  are  on 
different  planes  decay  exponentially.  Take  for  example 


GAa(R)  =  ^ 


J-lb  _  _ Jo(A.r)  X  U2b _ 

2ji  j[llbl3U3b  +  Uib  COtKbibUib)]  U2b  COSh{u2b{b2b  -  bib)) 

I  +  [libi2uib  +  |ib23UibU3b  coth(bibUib)]  Sinh(u2b{b2b  -  bib)) 


(4-16) 


which  is  the  Green’s  function  for  the  magnetic  vector  potential  on  interface  2b  for  a  point 
source  on  interface  3b  with  a  radial  separation  R.  Again,  as  A,  -4  «>  the  integrand 


becomes 


^  Itbt3  A.  +  X  ]  X  exp{Xtb)  +  .|ibi2X^  +  ltb23X^.  exp(Xtb))  (4  - 17) 


which  further  reduces  to 


3s(XR-|) 


It  I  exp(-Xtb) 

4^  jibl3  +  1  +  ^bl2  |ib23 


(4-18) 


and  the  dominant  factor  as  X  — >  ««>  is  exp[-Xtjjl.  Similar  results  are  found  for  all  other 
Green’s  functions  where  the  source  and  observer  are  on  different  planes. 


4.1.4  Behavior  as  R  —¥0.  To  determine  the  currents  using  the  moments  method 


outlined  in  chapter  ID,  it  is  necessary  to  perform  a  surface  integration  over  the  source 
(basis)  and  observer  (test)  areas.  When  the  source  and  observer  points  are  on  separate 
planes,  the  Green’s  functions  can  be  solved  for  a  discrete  series  of  R  values  between  R 
=  0  and  the  maximum  radial  separation  expected  in  the  antenna  design.  The  surface  inte¬ 
grations  are  then  accomplished  by  interpolating  the  values  of  the  various  Green’s  func¬ 
tions  for  different  /?’s  between  the  calculated  values  [51.  But  when  /?  =  0,  this  procedure 
fails  when  the  source  and  observer  are  on  the  same  plane. 

For  example,  when  R  =  0,  the  Bessel  function  in  (4  - 13)  becomes  1  and  (4  - 15) 
changes  to 

_ i  ^b23 _ 

^bl3  +  1  +  M-bl2  +  Hb23  (4  -  19) 
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as  X  — >  <».  There  is  no  decay  of  the  integrand  and  the  integral  blows-up,  as  expected, 
when  the  source  and  observer  are  at  the  same  point.  This  indicates  that  the  magnetic  vec¬ 
tor  potential  must  go  to  infinity  when  the  source  and  observer  are  at  the  same  point,  but 
the  source  current  is  still  finite. 

As  /?  — >  0,  a  different  approach  is  taken  to  apply  the  moments  method  to  Green’s 
functions  with  co-planer  source  and  observer  points.  An  asymptotic  approximation  to 
(4-  13)  is  made  by  first  splitting  the  interval  of  integration  at  some  and  doing  the 
change  of  variables  x  =  XR  to  obtain 


27C 


^  U2b  cosh(u2b  (b2b  -  bib))  +  )ib23  U3b  sinh(u2b  (b2b  -  bib))  ^ 

D^X) 


l^lb  I  Jo(x) _ 1  +  lib23 

L  ^  +  1  +  |ibl2  + 

jKXc 


db23 


■dx 


(4-20) 


where  X  is  chosen  large  enough  to  satisfy 


U2b  cosh(u2b  (b2b  -  bib))  +  M-b23  U3b  sinh(u2b  (b2b  -  bib)) 

Dag 

_ _ 1  +  lib23 _ 

^bl3  +  1  +  Ubl2  +  ^b23  (4-21) 

to  some  arbitrary  precision  (remember  U2j,  and  u^j^  are  also  functions  of  X^,).  From 
tables  [18]  it  is  known  that 


f 


Jo(x)dx  =  1 


(4  -  22) 


which  can  also  be  written  as 


I  Jo(x)dx  =  1  -  f  . 
JrXc  Jo 


Jo(x)dx 


(4  -  23) 


Using  (4  -  23)  in  place  of  the  second  term  in  (4  -  20)  then  yields 
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^A22(^)  *  I  Jq^X.R.)  a*  COSh(u2b  (b2b  ~  ^Ib))  M-b23  ^3b  sinh(u2b  (b2b  ~  ^Ib)) 


1  +  }^b23 _  f 

+  1  +  M-bl2  +  M-b23 


D^X) 

Jo(x)  ,  1  1  +  ^b23 

R  R  |Ibi3  +  1  +  |Ibl2  +  ^b23 


(4-24) 


The  first  and  second  integrals  in  (4  -  24)  are  now  over  finite  intervals  and  the  solutions 
finite  for  any  value  of  R.  If  the  Bessel  function  in  the  second  integration  is  expanded  into 
its  series  expression  and  integrated  term  by  term,  it  is  seen  the  value  of  the  integral  goes 
to  Xj,  as  R  — >  0.  The  surface  integration  of  the  last  term  in  (4  -  24)  is  accomplished  ana¬ 
lytically  for  the  moments  method  analysis.  Note  the  last  term  in  (4  -  24)  corresponds  to 
the  homogeneous  Green’s  function  [13]  for  an  unbounded  medium  of  permeability 


2(1  -*-}^b23) 

Hbl3  +  1  +  llbl2  +  I^b23 


This  asymptotic  approach  is  used  for  small  R  for  all  Green’s  functions  where  the 
source  and  observer  points  are  on  the  same  plane.  All  asymptotic  approximations  are 
given  along  with  the  regular  functions  in  Appendix  D.  The  asymptotic  method  is  not 
needed  for  Green’s  functions  with  source  and  observer  points  on  separate  planes,  since 
the  exponential  decay  of  these  integrands  is  independent  of  /?,  therefore,  the  integrals 
always  converge.  The  characteristics  of  the  integrands  for  all  the  Green’s  functions  are 
summarized  in  Table  4-1. 


4.2  Numerical  Evaluation  of  the  Green’s  Functions. 

With  the  characteristics  of  the  integrands  now  defined,  the  numerical  integration 
of  the  Green’s  functions  are  broken  up  into  two  general  cases  (see  Figure  4  -  4).  By 
breaking  up  the  integration  intervals,  the  unique  difficulties  are  isolated  and  handled 
separately.  The  imaginary  part  of  the  integral  in  both  cases  is  readily  calculated  using  the 
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Table  4-1  Green’s  Functions  Characteristics 


Behavior  as 

Integrand’s  rate 

Limit  of  integrand  as 

R^O 

X,  00  for  R  =  0 

'JA22 

Infinite 

^-m 

1 

1  +  llal2 

1 

1  +  £312 

G“22 

Infinite 

5,-172 

^bxx 

Ga22 

5,-172 

1  +  M-b23 

Infinite 

1  +  4bl2  +  ^b23  +  ltbl3 

k 

5,-172 

1  +  eb23 

Gq22 

Infinite 

1  +  ebl2  +  Eb23  +  ebl3 

/-.bxx 

Ga32 

Finite 

exp[-X.  tjj] 

NA 

Gq32 

Finite 

exp[-X.  tjj] 

NA 

^bxx 

Ga33 

5,-172 

1  +  ltbl2 

Infinite 

1  +  Itbl2  +  Itb23  +  4bl3 

h 

5,-172 

1  +  ebl2 

Gq33 

Infinite 

1  +  ebl2  +  eb23  +  £bl3 

G^^5 

Finite 

exp[-X  t{j] 

NA 

Gq23 

Finite 

exp[-X  tjj] 

NA 

x^axx 
;  ^Fll 

Infinite 

5,-172 

1 

■  ,-ia 

^mll 

Infinite 

5,-172 

1 

Gpii 

Infinite 

5,-172 

1 

G^ii 

Infinite 

5,-172 

1 

GI21 

Finite 

exp[-A.  bj^] 

NA 

G^i 

Finite 

expI-Xb^j] 

NA 

Gki 

Finite 

exp[-Xbu,] 

NA 

Ghi2 

Finite 

exp[-Xbj3] 

NA 

Gh12 

Finite 

expI-Xbjf,] 

NA 

JH13 

i 

Finite 

exp[-A.  bj,j] 

NA 

-60- 


IV.  Application  and  Results 


(b) 


Figure  4-4  General  integration  characteristics  for  the  Green's  functions. 


quadrature  routine  DQDAGS  from  the  IMSL  libraries.  The  real  parts  of  the  integrands 
are  piece-wise  continuous  due  to  the  integrable  singularity  at  X,  =  k^.  The  DQDAGP  rou- 
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tine  from  the  IMSL  libraries  is  used  successfully  to  integrate  over  region  I  in  both  cases 
using  X  =  ICq  as  the  location  of  the  singularity  [14].  Region  III  in  case  2  presents  no 
difficulties  and  DQDAGS  is  used  there.  Region  II  in  case  1  and  region  IV  in  case  2  rep¬ 
resent  the  part  of  the  integration  over  the  semi-infinite  interval,  and  the  Method  of  Aver¬ 
ages  [12]  described  below  is  used  in  this  region.  The  integration  intervals  for  these  re¬ 
gions  in  the  asymptotic  cases,  though  finite,  are  large  and  present  their  own  problems. 
The  final  difficulty  lies  in  region  H  for  case  2.  Some  integrands  contain  a  simple  pole 
and  the  integral  exists  only  in  the  Cauchy  principle  value  sense.  Once  the  location  of  the 
pole  is  found  using  methods  previously  described,  a  simple  routine  is  used  to  find  the 
principle  value  of  the  integral  between  ±A  as  shown  in  the  inset  of  Figure  4  -  4b. 

4.2.1  Taking  the  Cauchy  principle  value.  In  integrands  with  a  simple  pole,  the 
subroutine  cauchy  (see  Appendix  F)  is  used  to  find  the  principle  value  of  the  integral 
over  a  symmetrical  interval  about  the  the  pole.  Although  more  elegant  methods  exist  [12] 
to  find  the  principle  value,  the  method  used  here  simply  integrates  inward  to  the  pole 
starting  at  ±A  over  successively  smaller  intervals  (see  Figure  4  -  5)  until  the  solution  con¬ 
verges  to  some  defined  accuracy.  Since  the  pole  is  usually  close  to  (or  k2g),  the  best 
approach  in  determining  A  is  to  make  A  some  fraction  of  the  interval  between  k^jj  and 
the  pole  location.  The  expression 


Ab  = 


X.pb  -  k3b 
10 


(4-25) 


was  found  to  work  well  in  region  b. 

The  cauchy  subroutine  was  tested  (see  test  program  1,  Appendix  F)  on  the 
function  [12] 


■"f 


^^^^dX=  1.450590 
X- 1 


(4-26) 


using  A  =  0.1  and  returned  I  =  1.45058958.  As  long  as  the  exact  position  of  the  pole  is 
known,  this  procedure  converges  to  the  correct  solution  with  as  much  precision  as  the 
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Figure  4-5  Successive  integration  intervals  lor  Cauchy  principle  value. 


computer  is  capable  of. 

If  the  position  of  the  pole  is  approximated,  the  accuracy  of  the  solution 
deteriorates.  The  solution  converges  until  the  integration  interval  on  one  side  of  the  ap¬ 
proximated  pole  gets  too  close  to  the  actual  pole  position.  At  this  point  the  solution  starts 
to  diverge.  Table  4-2  gives  the  solutions  the  test  program  converged  to  for  several  ap¬ 
proximated  pole  positions.  Even  though  the  first  two  pole  positions  are  accurate  to  six 
significant  figures,  the  solutions  are  only  accurate  to  three.  It  took  pole  positions  accu¬ 
rate  to  twelve  figures  to  obtain  solutions  accurate  to  six  significant  figtires. 


Table  4-2  Solutions  With  Approximate  Pole  Positions 


Approximated  pole 

Solmion 

1.000001 

1.44963661 

0.999999 

1.45124023 

1.000000000001 

1.45058621 

0.999999999999 

1.45059361 
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The  solution  is  guaranteed  to  diverge  when  attempting  to  find  a  cauchy  principle 
value  by  integrating  over  successive  intervals  in  towards  an  approximate  pole.  To  insure 
convergence,  the  subroutine  cauchy  uses  what  may  be  called  a  reset  method.  The  size  of 
each  integration  interval  on  each  side  of  the  pole  in  Figure  4  -  5  is  half  the  previous  one 
(i.e.  Aj^  =  A/n  ).  At  the  integration  interval  the  solution  starts  to  diverge,  the  subroutine 
resets  the  solution  to  that  of  the  last  convergent  interval  and  continues  to  integrate  in  to¬ 
wards  the  pole  (see  Figure  4  -  6).  Except,  the  interval  is  again  halved  and  the  integrations 
get  no  closer  to  the  approximated  pole  than  the  point  where  the  solution  starts  to  diverge. 
The  subroutine  continues  to  halve  the  integration  intervals  until  the  desired  relative 
accuracy  is  reached.  Unfortunately,  as  shown  in  Table  4  -  2  the  answer  this  method  con¬ 
verges  to  is  not  very  accurate  unless  the  approximate  pole  position  is  much  more  accu¬ 
rate. 


interval.  Contribution  of  this  Inter¬ 
val  is  thrown  out,  interval  halved 
again,  but  started  at  same  point 

Figure  4*6  Reset  method  for  taking  Cauchy  principle  value. 


42.1.1  Calculating  the  residue.  If  an  integral  must  be  evaluated  in  the 
Cauchy  principle  value  sense,  it  also  has  a  residue  that  must  be  calculated.  Since  the  pole 
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location  is  approximated  and  cannot  be  factored  out  of  the  denominator  of  the  integrands, 
function  (4  -  4)  cannot  be  used  directly  to  find  the  residue.  If  the  part  of  the  denominator 
containing  the  pole  (D^  or  for  only  one  surface  wave  mode)  is  expanded  into  a  Tay¬ 
lor  series  about  the  pole,  the  residue  is  found  by 


residue 


=  lim  = 

^  ^  X,  dM  dtXp) 


(4-27) 


where  ng(Xp)  is  the  non-singular  part  of  the  integrand  and  d'(Xp)  is  the  first  derivative  of 
the  singular  part  of  the  denominator  evaluated  at  the  pole. 

For  example,  the  solution  to  Gq32{R)  for  any  value  of  R  is 


Jo(xr)x 


U2b 


Jo 


j  Jo(^b^)  ^pb  ti2b 


[u3b  +  Ubn  uib  tanh(bib  uib)] 

*  U2b  cosh(u2b  (b2b  -  bib)) 

ltbl3  ebl2  uib  +  ()tb23  *  P-bl3  Ebu)  uib 

,  .  +  dbi2  Uib  U3b  tanh(bib  Uib) 

\*  sinh(u2b  (b2b  -  bib)) 

|[u3b  +  ltbi3  uib  tanh(bib  Uib)]  1 

*U2bCosh(u2b(b2b-bib))  1 

Hbl3  £bl2  uib  +  {^b23  '  HbI3  ^bll)  U^b 
+  Kbi2  Uib  U3b  tanh(bib  Uib) 

\*  sinh(u2b  (b2b  -  bib)) 


dX 


(4-28) 


where 


T^P(-\  )_'!  ! .  (ebi2uib  +  eb23  uibU3btanh(bibUib))  , /U3b  .  U2b\„ 

D  ml  V/  -  '^pbj[tb - ^ ^ 

+  bib  U2b  sech^bib  uib)  +  tanh(bib  uib)  cosh(t  b  U2b) 

+  [tb(ebi3  U3b  +  Uib  tanh(bib  uib))  +  2  ebi2 

+  bib  eb23  U3b  sech^bib  uib)  +  eb23  tanh(bib  uib)  sinh(t  b  U2b)j 

(4-29) 

Note  u  jjj,  U2tj,  and  u^jj  are  all  evaluated  at  Xpj,  in  (4  -  29)  and  the  second  part  of  (4  -  28). 
Similar  expressions  arc  used  in  region  a  for  integrations  with  singularities  with 

D  m(^pa)  =  ^pa  ^  +  bi »  sech^bi  a  Ui  a)  + 


LU2 


ui, 


(4  -  30) 
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where  and  \X2^  are  evaluated  at 

42.2  Method  of  Averages.  The  Method  of  Averages  is  a  technique  developed  by 
Mosig  and  Gardiol  [12]  specifically  to  solve  Sommerfeld  integrals,  but  can  be  applied  to 
almost  any  integral  of  an  oscillating  function  from  a  finite  point  to  infinity.  Consider  the 
integral 


1  = 


Xpf(x)dX 


(4-31) 


where  f(X)  is  a  continuous  function  with  asymptotic  behavior 


lim  f(x)  =  C  X 

(4 . 32) 


where  a  =  -1/2  for  the  functions  in  this  paper  and  C  is  some  constant.  The  cases  with  ex¬ 
ponential  behavior  will  be  treated  later.  Integrals  like  (4-31)  can  be  evaluated  as  an 
oscillating  series  like  (4  - 12).  A  partial  series  Im  (m  =  1,  2,  .  .  .,  M)  for  (4  -  31)  is 
defined 


m  =  1,  2,  .  .  .,  M 


(4-33) 


where  X^^  are  the  successive  zeros  of  the  oscillating  function  greater  than  the  lower  limit 
of  integration  a.  Like  an  alternating  series,  the  difference  between  the  real  value  /  of  the 
integral  and  the  panial  series  value  li,  is  always  less  than  the  value  of  the  first  term  of  the 
series  neglected  [16].  Approximating  f(X)  with  the  first  term  of  its  Taylor  series  expan¬ 
sion  about  Xjjj,  this  first  neglected  term  can  be  integrated  analytically  and  the  prior  state¬ 
ment  written 

I-Im< 


2 

p 


fesl 

P 


(4  -  34) 


The  convergence  is  thus  determined  by  f(Xj^). 

A  new  sequence  is  then  defined  Im  (m  =  1,  2, . . .,  M-1)  by  taking  the  average  of 
two  successive  values  of  li,,  following  the  general  expression  [12] 
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m  =  1,...,  M-1,  n  =  1 . M-1 


(4-35) 


Taking  the  average  of  the  last  term  in  the  partial  series  li,  and  the  first  neglected  term,  the 
inequality  is  now  written 


I 


where 


I-Im<  ^ftxj 


^  -  iK  -  7t/p) 

Jt/p 


(4-36) 


(4-37) 


(Recall  sin(Xj^p  -  :c)  =  -  sin(3,j^p) ).  The  convergence  of  is  now  dependent  on  the  fu'st 
derivative  of  f(A,).  Subsequent  use  of  the  average  relationship  (4  -  35)  produces  new  se¬ 
quences  Im  that  converge  even  faster.  The  last  sequence  reduces  to  the  single  value 
which  is  closer  to  the  real  value  than  Im,  even  though  no  new  evaluations  of  the  integrand 


were  made  [12].  The  final  value  is  expressed  directly  in  terms  of  the  starting  sequence 
Imby 


M 


m=l 


(4  -  38) 


The  method  of  averages  can  also  be  applied  to  Bessel  functions  J^(>.p).  Using  the 
asymptotic  expression  for  the  Bessel  functions,  (4  -  34)  -  (4  -  36)  remain  valid  by  defin¬ 
ing  the  values  as  the  zeros  of  cos(Xp  -  Jt/4  -vk/2)  and  replacing  f(X)  with 

f(X)(2/7tXp)^/2  JJ2] 

It  should  be  noted  that  (4  -  35)  suggests  an  optimal  sequence  can  be  obtained 
by  taking  as  the  mean  points  between  the  zeros  of  cos(Xp).  The  sequence  li,  then 
converges  as  ffX)  and  one  less  averaging  operation  is  required  [12]. 

The  convergence  can  be  speeded  up  even  more  by  introducing  a  weighted  aver¬ 
age.  The  arithmetical  mean  (4  -  35)  is  replaced  by  a  weighted  average  where  more 
weight  is  given  to  the  values  of  li,  closest  to  I  [12].  This  is  expressed  as 
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jn-^1  ^  V'm  Im  +  Wm-H  Im^l  n=l,...,  M-1 

+  (4-39) 


For  n  =  1 ,  the  optimal  value  of  weights  is  found  from  (4  -  34)  to  be 

\« 


Wm  = 


h 

\lj 


(4  -  40) 


where  a  is  the  asymptotic  exponent  (decay)  of  the  function  f(A.).  After  each  successive 
averaging  operation,  the  order  of  the  function  controlling  convergence  decreases  by  one 
unit  [12],  so  (4  -  40)  is  generalized  to 

/.  \a  +  1  -  n 


,n  _  I  ''■1 


Wm  = 


(4-41) 


The  optimal  series  for  Bessel  functions  Jv(R^)  is 


(4-42) 


When  applied  to  an  actual  problem,  the  Method  of  Averages  is  applied  iteratively 
until  a  specified  relative  accuracy  is  achieved.  The  sequence  Im  is  built  up  one  element  at 
a  time  by  integrating  over  successive  intervals  and  the  best  solution  found  and  com¬ 
pared  to  the  previous  best  solution.  In  this  manner,  only  the  last  m  value  of  each  se¬ 
quence  C  plus  the  value  of  the  next  integration  from  to  are  needed  to  calculate 
the  next  best  solution  \  An  example  of  this  process  is  shown  in  Figure  4  -  7.  The 
circled  values  are  the  only  ones  needed  for  the  next  iteration.  The  next  iteration  is  started 
by  integrating  the  function  over  the  next  interval  to  obtain 

l6  =  l5+[  gWdX 

Aj 


(4  -  43) 


where  g(X)  is  some  decaying,  oscillating  function.  By  applying  weighted  averaging  to 
this  value  with  the  previous  results,  a  new  solution  is  found  that  is  more  accurate  than  the 
previous  one.  This  cycle  is  repeated  until  the  difference  between  the  new  and  previous 
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l6 


Integration  over  new 
interval  +I5 


solutions  is  less  than  some  desired  relative  accuracy. 
The  subroutine  bavg  was  tested  on  the  function 


(4-44) 


using  the  optimal  sequence  =  (m  +  l/4)ji  (see  Test  Program  2  in  Appendix  F).  Each 
half-period  was  evaluated  using  the  DQDAG  subroutine  from  the  IMSL  library  which 
uses  an  adaptive  Gauss-Kronrod  rule  with  10-21  point  quadrature  to  estimate  the  inte¬ 
gral  [14].  The  results  are  presented  in  Figure  4-8.  Curve  A  represents  the  sequence  1^, 
(no  averaging)  and  shows  the  slow  convergence  expected  of  a  function  that  decays  as 
Curve  B  depicts  the  accelerated  convergence  after  straight  averaging  was  used  to 
obtain  l"  (n  =  m).  Curve  C  was  produced  using  weighted  averaging.  The  fastest 
conve.  ^.jnce  is  obtained  by  this  case,  reaching  an  accuracy  of  10'^  or  better  after  only  6 
evaluations. 
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Figure  4  -  8  Relative  error  in  numerical  evaluation  of  (4  -  44)  (E  =  |I*- 1 1/ 1  where  I*  is  the 

numerical  solution). 

4.22.1  Using  the  Method  of  Averages  with  exponentially  decaying  func¬ 
tions.  Mosig  and  Gardiol  focused  on  a  single  layer  antenna  where  all  the  source  and  ob¬ 
server  points  are  on  a  single  interface  [12].  Hence,  all  their  Green’s  functions  decayed  as 
^-1/2  Method  of  Averages  they  developed  addressed  only  the  case  of  algebraic 
decay  and  did  not  consider  the  case  of  exponential  decay  obtained  when  the  source  and 
observer  are  on  separate  planes. 

The  Method  of  Averages  is  extended  here  to  exponentially  decaying  functions 
using  similar  reasoning  as  above.  First,  consider  the  integral  (4  -  31),  where  f(^)  is  a 
continuous  function  with  asymptotic  behavior 

lim  f{x)  =  C  exp(-px) 

X-K«>  (4  -  45) 

As  before,  a  partial  series  li,  is  formed  by  using  (4  -  33)  that  is  as  accurate  as  the  magni¬ 
tude  of  the  first  term  neglected.  Instead  of  using  the  first  term  of  the  Taylor  series  of  f(X) 
to  approximate  the  neglected  term,  the  asymptotic  approximation  of  (4  -  45)  is  used  to 
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obtain 

I  -  Im  <  exp(-pXj[exp(-p7r/p)  +  l] 

P  +  (4  -  46) 

and  the  rate  of  convergence  is  determined  by  the  exponential  decay  as  expected. 

'y  1 

Calculating  the  new  sequence  Im  by  averaging  successive  values  of  Im  as  before, 
the  inequality  is  now  written 

I  -  Im  <  exp(-pXm)  sinh(pJi/p) 

P  +  P^  (4  -  47) 

2  1 

This  result  shows  Im  may  or  may  not  converge  faster  than  Im,  depending  on  the  value  of 
p7t/p.  Therefore,  straight  averaging  does  not  guarantee  the  fastest  convergence. 

This  dilemma  is  resolved  by  using  a  weighted  average.  The  weight  function  is 
found  similarly  as  before 

Wm  =  =  exp[p(^m  -  >-i)] 

exp(-pXmj  (4  -  48) 

Note,  the  weight  function  does  not  depend  on  the  level  of  averaging  as  in  (4  -  41),  only 

the  value’s  position  in  the  partial  series.  If  (4  -  48)  is  used  to  calculate  a  weighted  aver- 

age  series  Im,  the  convergence  is  found  to  be 

I-Im<0  (4-49) 

Of  course,  the  asymptotic  expression  for  f(X)  is  not  perfect  and  the  right-hand-side  of 

(4  -  49)  will  not  be  exactly  zero.  What  this  does  show  is  that  weighted  averaging 

converges  faster  than  the  partial  series  and  straight  averaging  independent  of  pJt/p. 

The  bavg  subroutine  was  tested  on  the  exponentially  decaying  function  [18] 

1=1  ex  p(-pX)jo^pX)dX 

Jo  (4  -  50) 

where  I  =  l/sqrt(2)  for  P  =  p  =  1  and  I  =  l/sqn(1.01)  for  P  =  0.1,  p  =  1  (see  Test  Program 
3  in  Appendix  F).  In  both  Figure  4  -  9a  and  Figure  4  -  9b,  curve  A  represents  the  partial 
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Figure  4  -  9  Method  of  Averaging  applied  to  exponentially  decaying  function  (4  -  50)  for 

(a)p  =  p=l:(b)p-0.l.p.i. 


series  Im,  curve  B  represents  the  results  for  straight  averaging,  and  curve  C  corresponds 
to  the  results  for  weighted  averaging.  Note  in  both  cases  weighted  averaging  produces 
the  best  convergence. 

4. 2. 2. 2  Integration  of  exponentially  decaying  function  as  R  —>0.  In  theo¬ 
ry,  as  long  as  ^  >  0  the  integrand  oscillates  due  to  the  Bessel  function  and  an  alternating 
series  that  converges  is  formed.  In  reality,  as  0  the  periods  of  the  Bessel  function 
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become  very  long,  making  numerical  integration  difficult  over  the  half  periods  used  in 
the  Method  of  Averages.  Attempting  to  do  so  typically  results  in  overflow  errors  in  the 
integration  subroutines.  Since  the  integrand  still  decays  exponentially  independent  of  R, 
it  is  easier  to  simply  integrate  over  a  series  of  finite  intervals  and  sum  the  result  until 
some  desired  relative  accuracy  is  satisfied. 

4.2.3  Integration  for  the  asymptotic  case.  When  /?—>  0  for  the  Green’s  func¬ 
tions  where  the  source  and  observer  are  on  the  same  plane,  the  integral  blows  up.  This 
behavior  and  a  technique  to  deal  with  it  were  described  in  section  4.1.4.  To  implement 
this  technique,  two  problems  must  be  surmounted.  First,  the  critical  value  X.  must  be 
found  where  the  integrand  can  be  approximated  by  its  asymptotic  form.  Second,  a  meth¬ 
od  to  integrate  the  original  integrand  and  the  function  Jq(x)  from  0  to  X^  and  0  to  RX.^,  re¬ 
spectively,  must  be  developed. 

4.2.3. 1  Finding  X^.  The  critical  value  is  found  by  •'pplying  a  bisection 
method  [19].  A  function  of  the  integrand  with  R  =  0  (Bessel  function  is  1)  minus  the  as¬ 
ymptotic  limit  of  the  integrand  is  used  to  find  the  X^  where  this  function  equals  some 
stopping  tolerance. 

For  example,  the  function  used  to  find  X„  for  Ga22  is 

V 

^  U2b  cosh(u2b  {b2b '  bib))  +  lib23  U3b  sinh(u2b  (b2b  -  bib)) 

A,r - - - 

_ ^  l^b23 _ 

M-bl3  +  I  +  llbl2  +  ltb23  (4-51) 

As  X^  — >  o®,  (4-51)  0.  The  subroutine  brackets  the  point  X^  where  the  function 

equals  the  stopping  tolerance  (10'^  was  found  to  work  well),  defining  points  where  the 
function  is  higher  and  lower  than  the  stopping  value.  With  the  bracketing  points  estab¬ 
lished,  it  is  known  the  value  X^  exists  somewhere  in  the  bracketed  interval.  The  interval 
is  then  bisected  repeatedly,  keeping  the  point  where  the  function  equals  the  stopping  cri¬ 
teria  in  the  present  interval.  Since  the  stopping  criteria  is  given  and  the  size  of  the 
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original  interval  is  defined  in  the  subroutine,  the  number  of  iterations  to  make  using  the 


bisection  method  is  known  by 

n  =  logj— 
A 


(4-52) 


where  Aq  is  the  size  of  the  original  interval  and  A  is  the  given  stopping  tolerance  [19]. 

Although  any  value  greater  than  the  found  with  this  method  is  sufficient  to  de¬ 
fine  a  point  where  the  asymptotic  approximation  can  be  made,  it  is  important  not  to  make 
any  larger  than  necessary  to  achieve  the  desired  accuracy.  Remember,  the  intervals  of 
two  numerical  integrations  depend  upon  the  value  of  X^;  the  larger  these  intervals  are,  the 
more  difficult  and  time  consuming  the  integrations  become. 

4.2.3. 2  Integrating  to  and  RX^.  To  complete  the  analysis  for  the  as¬ 
ymptotic  case,  the  original  integrand  and  Jq(x)  must  be  numerically  integrated  from  0  to 
X^  and  0  to  RX^,  respectively  (see  (4  -  24)  as  an  example).  The  difficulties  with  possible 
singularities  and  the  complex  behavior  of  the  original  integrand  for  X  <  k^jj  (or  X  <  k2a) 
are  handled  as  shown  in  the  previous  sections  of  this  chapter.  Since  X^  is  usually  very 
large,  both  the  original  integrand  and  Jq(x)  will  oscillate  many  times  over  their  respective 
intervals  for  even  a  very  small  value  of  R.  Applying  normal  quadrature  methods  to  these 
functions  over  their  entire  integration  intervals  would  be  difficult.  Instead,  the  intervals 
are  divided  into  a  number  of  finite  sized  sub-intervals  no  greater  than  one  half  period 
long.  The  functions  are  then  fairly  smooth  over  each  sub-interval  and  normal  quadrature 
subroutines  can  be  applied  over  the  sub-intervals  and  the  results  added  together. 


4.3  Sample  Results  for  the  Green’s  Functions. 

The  numerical  techniques  developed  above  were  applied  to  several  of  the  Green’s 
functions  calculated  in  chapter  ni.  Results  were  obtained  for  the  functions  Gam.  ^’q22, 
^q23,  and  Gliri  all  the  Green’s  function  necessary  to  describe  the  tangential  elec¬ 
tric  fields  at  interface  2b  for  all  the  sources  in  region  b.  The  functions  G^3,  Gq23,  and 
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Ge2i  were  solved  for  25  values  of  R  between  0  and  0.05  m.  The  functions  Ga22  and 
Gq22  were  solved  for  25  points  between  0.001  and  0.05  m  and  their  asymptotic  functions 
were  solved  for  25  point  between  0  and  0.005  m.  A  quadratic  distribution  [12]  of  points 
was  used  so  a  higher  concentration  of  data  points  would  be  taken  near  R^jj^  where  the 
functions  change  more  rapidly.  The  quadratic  distribution  may  be  described  by 


(n^-  l) 


(4  -  53) 


where  N  is  the  total  number  of  points  to  evaluate.  All  functions  were  evaluated  for  the 
following  parameters; 


frequency  =  4  Ghz 
b^jj  =  0.0016  m 
^2b  "  0.0048  m 
eib  =  2.5eo 

^2b  =  5  ^0 

^3b  =  ^0 

^^Ib  =  ^^2b  =  ^3b  =  ^^0 

For  these  values  the  pole  was  found  to  be  \b  =  88.73775940774  1/m.  The  critical  value 
for  G^2  was  =  88126.1171933  1/m  and  the  critical  value  for  Gq22  was  X-  = 
88465.924790  1/m.  Both  critical  values  were  found  for  a  stopping  tolerance  of  10'^. 

4.3.1  Results  for  Ga22  ond  Gq22.  Figure  4  -  10a  and  c  plot  the  real  and 
imaginary  parts  of  G^2-  The  real  part  goes  to  infinity  as  expected  and  the  imaginary 
part  is  finite  for  all  values  of  R.  Figure  4  -  10b  is  a  more  detailed  plot  of  the  real  values 
of  G^2  for  R  >  0.01  ra.  This  plot  shows  the  beginning  of  the  oscillatory  natiite  of  the 
Green’s  functions.  The  real  and  imaginary  parts  of  Gq22  are  plotted  in  Figure  4  -  1  la 
and  b,  respectively. 

Figure  4  -  12a  is  a  plot  of  the  real  part  of  the  finite  functions  (the  sum  of 
the  two  integrals  in  (4  -  24))  of  the  asymptotic  expression  for  G^2-  Note  the  results  are 
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Figure  4  - 10  Plots  of  Gas:  (a)  real  part,  (b)  oscillatory  nature  of  real  part,  ;c)  imaginary  part. 
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finite  for  ^  =  0.  Figure  4  -  12b  is  a  comparison  the  complete  asymptotic  expression 
(4  -  24)  and  the  normal  expression  (4-13)  of  G^2  over  the  points  R  where  the  evalua¬ 
tions  overlapped.  The  two  curves  are  practically  the  same.  The  noticeable  difference  is 
mainly  due  to  the  fact  that  data  was  obtained  over  more  points  for  the  asymptotic  case 
than  the  normal  case  for  this  interval.  Similar  results  for  Gq22  are  plotted  in 
Figure  4-13. 

4.3.2  Results  for  G^  and  Gq23.  The  Green’s  function  G^3  is  plotted  in 
Figure  4-14.  Although  the  real  part  of  the  function  becomes  very  large  as  /?  ->  0,  it  re- 
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Figure  4-12  (a)  Real  part  of  asymptotic  functions  for  Gaz.  (b)  comparison  of  results  for 
asymptotic  [A]  and  normal  [B]  representations  of 


mains  finite.  The  function  also  oscillates  like  G^2  as  R  increases,  though  this  is  not 
seen  in  the  plot  due  to  the  scale.  The  real  and  im  .aginary  parts  of  Gq23  are  plotted  in 
Figure  4-15. 

4.3.3  Vector  plots  of  Ge2i-  The  Green’s  function  Ge21  represents  the  tangential¬ 
ly  electric  field  at  interface  2b  for  a  magnetic  surface  current  source  on  the  ground  plane 
and  is  actually  a  dyadic.  The  results  from  three  different  Sommerfeld  integrals  (see  Ap¬ 
pendix  D)  are  used  to  calculate  the  vector  quantities  for  x-  and  y  -directed  magnetic  cur- 
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(a) 


b 

Figure  4-13  (a)  Real  part  of  asymptotic  functions  for  G^,  (b)  comparison  of  results  for 
asymptotic  [A]  and  normal  [B]  representations  of  G<^2. 


rents.  Numerical  and  graphical  results  for  these  individual  integrals  are  presented  in  Ap¬ 
pendix  F.  Assuming  an  infinitesimal,  unit  strength  (1  volt/m),  x-dirccted  magnetic  sur¬ 
face  current  located  at  the  origin;  plots  can  be  generated  depicting  the  relative  electric 
field  strength  and  direction  for  this  source  at  a  number  of  discrete  points.  Figure  4  -  16a 
represents  the  real  part  of  the  electric  field  on  interface  2b  at  0.01  m  intervals  between  the 
values  X  =  -0.05  m  to  x  =  0.05  m  and  y  =  -0.05  m  to  >  =  0.05  m  with  the  source  at  x  = 
y  =  0  on  the  ground  plane.  The  length  of  each  arrow  represents  the  relative  magnitude  of 
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Figure  4-14  Plots  of  Gas;  (a)  real  part,  (b)  imaginary  part. 


the  field  at  that  point  compared  to  the  greatest  magnitude  plotted.  Several  values  in  the 
center  of  Figure  4  -  16a  were  removed  because  their  magnitudes  were  so  high  the  relative 
sizes  of  the  other  arrows  were  too  small  to  reveal  any  detail.  Those  values  removed 
produced  very  large  arrows  pointing  to  the  top  of  the  figure.  Figure  4  -  16b  represents 
the  imaginary  part  of  the  electric  field. 

4.3.4  Interpolation  of  data.  The  calculation  of  Ge2i  at  the  specific  points 
for  Figure  4-16  required  the  interpolation  of  data  from  calculated  points.  The  integrals 
that  compose  Ge2i  were  solved  for  a  quadratic  distribution  of  25  values  of  R  between  0 
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b 

Figure  4*15  Plots  of  Gqzs;  (a)  real  part,  (b)  imaginary  part. 


and  0.05  m.  Obviously,  all  the  points  plotted  in  Figure  4  -  16  do  not  correspond  with  the 
values  of  R  where  these  integrals  were  solved.  It  would  be  too  time  consuming  to  solve 
these  integials  for  each  specific  value  of  R  necessary  for  these  calculations  (or  for  the 
moment  method  calculations),  therefore,  interpolation  was  used  to  find  values  for  the  in¬ 
tegrals  between  the  calculated  values. 

A  simple  polynomial  averaging  interpolation  [12]  was  used  to  calculate  the 
data  for  Figure  4-16.  The  method  fits  two  second  order  polynomials  to  the  four  points 
immediately  surrounding  the  interpolated  point  R  (see  Figure  4  -  17),  one  equation  for 
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Figure  4  - 16  Vector  plot  of  Gea  for  x-directed  source;  (a)  real,  (b)  imaginary  parts. 
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R-  Rj,  and  R-_|_j  and  the  other  for  Rj,  Rj+p  and  Rj^2-  polynomial  equations 

are  solved  for  R  and  the  results  averaged  to  obtain  the  interpolated  value.  This 
technique  was  executed  in  Mathematica  with  the  polyAvg.m  package  and  the  results 
plotted  in  Figure  4  -  16  using  the  vectorPIot.m  package,  both  listed  in  Appendbc  F. 
Polynomial  averaging  interpolation  is  also  well  suited  for  interpolating  data  to  use  in  the 
surface  integrations  for  the  moments  method. 


Figure  4-17  Polynomial  averaging. 


4.3.5  Determining  error  estimates.  Since  the  integrals  used  to  calculate  the 
Green’s  functions  must  be  solved  numerically,  it  is  important  to  have  some  idea  of  what 
relative  errors  are  inherent  in  the  solurions.  The  quadrature  subroutines  used  from  the 
IMSL  libraries  return  an  absolute  error  estimate,  where  the  difference  between  the  nu¬ 
merically  derived  result  and  the  actual  solution  for  the  integral  is  less  than  or  equal  to  the 
absolute  error  estimate  [14].  The  custom  integration  subroutines  written  for  this  research 
return  the  difference  between  the  results  of  the  last  two  iterations  as  the  absolute  error  es¬ 
timate.  The  absolute  error  estimates  obtained  for  the  integrations  over  the  different 
intervals  (see  Figure  4  -  4,  page  61)  are  then  summed  and  divided  by  the  final  result  to 
provide  an  overall  relative  error  estimate.  The  relative  error  estimate  reflects  the 
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precision,  or  significant  digits,  of  the  result.  For  example,  a  result  with  a  relative  error 
estimate  of  10'^  has  a  precision  of  6  significant  digits.  Unfonunately,  the  error  estimates 
for  the  principle  value  integrals  cannot  take  into  account  errors  generated  by  the  approxi¬ 
mated  pole  position.  For  this  reason,  it  is  imponant  that  the  accuracy  of  the  pole 
positions  be  much  greater  than  the  requested  relative  error  for  the  principle  value 
integrals. 
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A  theoretical  model  for  an  aperture  fed  stack-patch  microstrip  antenna  based  on  a 
full- wave  analysis  has  been  developed.  The  Green’s  functions  required  in  the  mixed 
potential  integral  equations  (MPIEs)  were  constructed  and  numerical  techniques  to 
evaluate  the  Green’s  functions  found.  A  moments  method  was  discussed,  but  not 
actually  implemented,  to  evaluate  the  MPIEs  for  the  electric  and  equivalent  magnetic 
currents  of  the  antenna.  The  remainder  of  this  chapter  answers  the  research  questions 
posed  in  the  introduction  of  this  thesis,  suggests  improvements  and  possible  applications 
for  this  research,  and  discusses  some  of  the  general  observations  made  during  the  study 
and  solution  process  involved. 

5.1  Answers  to  Research  Questions. 

5.1.1  What  are  the  Green’ s  functions  necessary  to  solve  for  the  various  currents 
of  the  antenna?  The  Green’s  functions  found  in  chapter  III  and  detailed  in  Appendix  D 
are  a  continuum  (integration)  of  elementary,  cylindrical  wave  functions,  generated  by  an 
infinitesimal  source  at  the  various  interfaces  of  a  stacked-patch  aperture  fed  microstrip 
antenna.  These  type  of  integrals  are  generally  referred  to  as  Sommerfeld  integrals.  The 
Green’s  functions  account  for  different  observation  frequencies,  dielectric  layer 
thicknesses,  dielectric  permeabilities  and  permittivities. 

5.1.2  What  are  the  mathematical  characteristics  of  the  integrands  in  the  Green’s 
functions?  The  integrands  of  the  Green’s  functions  are  decaying  oscillatory  functions 
with  semi-infinite  (0  to  <»)  intervals  of  integration.  All  integrands  are  complex  for  X  <  Icq 
and  real  for  X  >  Icq,  where  X  is  the  variable  of  integration  and  Icq  is  the  frce-space  wave 
number.  The  integrands  of  the  Green’s  functions  with  source  and  observer  points  on  sep- 
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arate  planes  decay  exponentially  and  the  integrals  remain  finite  as  /?  — >  0,  where  R  is 
the  radial  distance  between  the  source  and  observer  points.  The  integrands  of  the 
Green’s  functions  for  co-planer  source  and  observer  points  decay  algebraically  as 
and  the  integrals  blow-up  as  R  ^0.  An  alternate  asymptotic  expression  was  found  for 
small  R  for  these  cases. 

The  cut-off  frequencies  for  the  surface  waves  are  determined  by  the  zeros  of  the 
denominators  in  certain  integrands.  The  zeros  define  simple  poles  in  these  integrands 
and  require  these  integrals  be  evaluated  in  the  Cauchy  principle  value  sense.  The  zeros 
of  the  denominators  are  located  accurately  using  numerical  methods  and  the  principle 
values  of  these  integrals  evaluated.  With  a  proper  choice  of  material  parameters  and  op¬ 
erating  frequency,  only  one  pole  is  present  in  these  integrals  (only  one  surface  wave 
mode  propagates),  greatly  easing  the  evaluation  of  these  integrals. 

5.1.3  How  can  the  Green’s  functions  be  evaluated  numerically?  The  integrals  in 
the  Green’s  functions  are  evaluated  with  a  combination  of  commercially  available 
quadrature  subroutines  from  the  IMSL  libraries  and  custom  subroutines  written  to  handle 
the  specific  difficulties  of  these  integrals.  A  subroutine  is  written  to  find  the  principle 
value  of  an  integrand  with  a  simple  pole,  where  the  location  of  the  pole  is  approximated. 
A  subroutine  to  evaluate  an  integrand  from  some  finite  point  to  infinity  using  the  Method 
of  Averages  is  developed  and  the  original  technique  proposed  by  Mosig  and  Gardiol  [12] 
is  expanded  to  also  work  with  exponentially  decaying  integrands. 

As  /?  ->  0,  the  Method  of  Averages  brakes  down  and  alternate  methods  found. 
Those  integrals  that  are  finite  at  /?  =  0  (integrand  decays  exponentially)  are  integrated 
over  a  series  of  small,  finite  intervals  and  these  results  summed  until  a  desi’^ed  accuracy, 
determined  by  the  user,  is  met  Asymptotic  approximations  for  those  integrals  that  are 
infinite  at  ^  =  0  (integrand  decays  as  X'^^)  are  developed  and  special  subroutines  to 
integrate  the  resultant  oscillating  functions  over  large  intervals  ’.vritten. 
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5.1.4  What  methods  can  be  used  to  reduce  computational  time?  The  main  goal 
of  this  research  was  to  construct  the  Green’s  functions  and  demonstrate  the  resultant 
integrals  could  be  numerically  evaluated.  Computational  speed  was  a  secondary  concern, 
but  some  efficient  techniques  were  developed.  All  integrands  have  a  real  part  over  the 
interval  A.  =  0  to  ««  and  an  imaginary  part  from  A,  =  0  to  k^.  Numerical  integration  rou¬ 
tines  are  applied  separately  to  the  real  and  imaginary  parts  of  the  integrand  over  their 
respective  intervals.  Thus  saving  time  by  not  integrating  over  an  interval  where  the 
imaginary  part  of  the  integrand  is  zero. 

It  is  shown  the  integral  of  a  decaying,  oscillating  function  from  a  finite  point  to 
infinity  always  converges.  By  integrating  over  the  half-periods  of  the  oscillating  func¬ 
tion,  a  decreasing,  oscillating,  series  can  be  formed.  A  partial  sum  is  then  taken  that  is  as 
accurate  as  the  first  term  of  the  series  neglected.  The  rate  of  convergence  is  increased  by 
averaging  the  results  of  the  partial  series.  The  convergence  is  speeded  up  even  more  by 
applying  the  weighted  averaging  scheme  of  the  Method  of  Averages.  It  is  further  shown 
that  for  exponentially  decaying  integrands,  weighted  averaging  is  consistently  the  most 
efficient  method  to  use. 

The  use  of  the  asymptotic  approximations  should  be  restricted  to  cases  of  small 
R.  Although  the  asymptotic  approximations  can  be  used  for  any  value  of  R,  the  larger  R 
becomes  the  more  oscillations  appear  in  the  integration  intervals.  The  intervals  must 
then  be  divided  into  more  sub- intervals,  so  the  functions  are  “sufficiently  smooth”  for  the 
standard  quadrature  integration  subroutines  to  work;  thus  increasing  computation  time 
greatly. 

The  moments  method  analysis  requires  surface  integration  of  the  Green’s  func¬ 
tions  over  the  source  and  observer  coordinates,  p'  and  p  respectively.  The  integrands  of 
the  Green’s  functions  are  dependent  on  the  radial  separation  of  the  source  and  observer 
points  ( /?  =  Ip'  -  pi ).  To  numerically  evaluate  the  integrals  of  each  Green’s  function  at 
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every  value  of  R  required  in  these  surface  integrations  would  be  extremely  time  con¬ 
suming.  It  is  more  efficient  to  evaluate  the  integrals  in  the  Green’s  functions  over  a  finite 
number  of  R  values  and  then  interpolate  the  values  of  the  integrals  between  these  solved 
values.  It  was  shown  polynomial  averaging  may  be  an  acceptable  interpolation  tech¬ 
nique.  The  surface  integrations  can  then  be  done  using  the  interpolated  values  of  the 
Green’s  functions. 

5.1.3  How  can  the  solutions  of  the  Green’s  functions  be  used  in  a  moments 
method  solution  for  the  currents  of  the  antenna?  When  the  source  and  observer  are  in 
the  same  plane,  both  non-asymptotic  and  asymptotic  expressions  for  the  Green’s  func¬ 
tions  must  be  used.  The  surface  of  each  conductor  or  aperture  must  be  divided  into  a 
number  of  rectangular,  overlapping  current  cells  for  the  moments  method.  Surface  inte¬ 
grations  of  the  product  of  the  Green’s  functions  and  test  or  basis  functions  over  the 
surfaces  of  the  current  cells  can  then  be  accomplished.  When  the  source  and  observer 
points  are  in  separate  cells,  the  radial  separation  R  is  large  enough  to  use  the  interpolat¬ 
ed  values  from  non-asymptotic  expressions  of  the  Green’s  functions.  When  the  source 
and  observer  points  are  in  the  same  current  cell,  the  interpolated  values  from  the  asymp¬ 
totic  expressions  should  be  used  along  with  the  analytic  evaluation  of  the  1/R  term.  The 
dimensions  of  the  individual  current  cells  should  be  small  enough  to  make  the  evaluation 
of  the  integrals  in  the  asymptotic  expression  reasonably  efficient.  Using  the  non- 
asymptotic  expression  when  the  source  and  observer  are  in  separate  cells,  and  the 
asymptotic  expression  when  they  are  in  the  same  cell,  avoids  switching  between  the  two 
expressions  of  the  Green’s  functions  in  the  middle  of  a  surface  integration.  When  the 
source  and  observer  points  are  on  separate  planes,  no  asymptotic  approximations  are 
necessary  and  the  surface  integrations  can  be  accomplished  using  the  interpolated  values 
of  the  Green’s  functions.  These  methods  can  then  be  used  to  fill  the  sub-matrices  and  so¬ 
lution  matrix  of  the  moments  method.  The  solution  for  the  currents  of  the  antenna  can 
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then  be  found  using  normal  matrix  solving  techniques. 

5.2  Future  Improvements  and  Applications. 

5.2.1  Improvements.  Much  remains  to  be  done  to  complete  this  line  of  research. 
Only  five  of  the  Green’s  functions  were  coded  and  evaluated,  the  remaining  functions 
must  be  written  as  Fortran  functions  and  the  necessary  subroutines  applied  to  evaluate 
them.  A  program  must  also  be  written  to  implement  the  moments  method  to  solve  for  the 
currents  of  the  antenna.  It  will  then  be  possible  to  compute  the  antenna’s  input  imped¬ 
ance,  resonant  frequency,  bandwidth,  and  other  characteristics  as  they  vary  with  changes 
in  dielectric  permeability,  permittivity,  and  thickness;  and  changes  in  aperture  and  patch 
sizes,  and  relative  positions. 

The  theoretical  model  must  be  validated  against  experimental  results.  Actual  mi¬ 
crostrip  antennas  should  be  constructed  and  tested,  and  the  results  compared  to  those  ob¬ 
tained  with  this  model.  If  experimental  testing  cannot  be  accomplished,  it  may  be  possi¬ 
ble  to  compare  results  of  this  theoretical  model  with  results  from  some  kind  of  finite-ele¬ 
ment  simulation. 

Funher  analysis  of  the  Green’s  functions  is  necessary.  The  Green’s  functions  de¬ 
veloped  in  this  paper  are  applied  only  to  the  tangential  fields  at  the  various  dielectric  and 
ground  plane  interfaces  of  the  antenna.  In  any  practical  antenna  analysis,  it  is  desirable 
to  know  what  far-field  pattern  the  antenna  generates.  To  accomplish  this,  it  will  be  nec¬ 
essary  to  find  asymptotic  forms  of  the  Green’s  functions  for  large  R  and  z.  Once  the 
currents  of  the  antenna  have  been  found  with  the  near-field  Green’s  functions  and  mo¬ 
ments  method,  the  currents  can  be  used  with  the  far-field  Green’s  functions  to  determine 
the  far-field  pattern  of  the  antenna. 

The  effect  of  surface  waves  should  be  funher  explored.  Surface  waves  will  effect 
both  the  far-field  of  the  antenna  and  the  radiation  efficiency.  Power  radiated  as  surface 
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waves  will  decrease  radiation  efficiency.  Certain  configurations  of  dielectric 
permittivities,  permeabilities,  or  thicknesses  may  decrease  the  power  radiated  in  surface 
waves  and  increase  radiation  efficiency,  or  enhance  the  far-field  pattern  of  the  antenna. 

The  numerical  integration  routines  can  possibly  be  improved.  The  quadrature  in¬ 
tegration  subroutines  from  the  IMSL  libraries  were  used  simply  for  convenience.  With 
further  study  of  the  integrands  in  the  Green’s  functions,  it  may  be  possible  to  develop 
more  efficient  integration  routines  or  find  more  efficient  ways  to  apply  the  IMSL 
subroutines. 

A  very  simple,  basic  approach  was  taken  to  find  the  Cauchy  principle  value  of 
those  integrals  with  simple  poles  and  better  methods  may  exist.  Mosig  and  Gardiol  [12] 
mentioned  two  possible  methods,  folding  around  the  pole  and  extraction  of  the  singulari¬ 
ty.  It  is  also  desirable  to  obtain  a  better  estimate  of  the  error  in  the  principle  value  due  to 
the  approximation  of  the  pole  location. 

5.2.2  Ideas  for  future  applications.  The  most  desirable  goal  of  this  line  of  re¬ 
search  is  to  develop  a  set  of  design  criteria  for  an  apenure  fed  stacked-patch  microstrip 
antenna.  Once  the  model  proposed  in  this  paper  is  complete,  it  should  be  possible  to 
analyze  the  effect  different  material  parameters  have  on  antenna  performance.  By  find¬ 
ing  out  how  such  characteristics  as  resonant  frequency,  bandwidth,  radiation  efficiency, 
radiation  resistance,  and  far-field  pattern  change  with  changes  in  different  antenna  pa¬ 
rameters,  it  should  be  possible  to  determine  what  parameters  produce  the  best  design  for 
a  given  application. 

The  results  of  this  research  are  by  no  means  confined  to  the  study  of  aperture  fed 
stacked-patch  microstrip  antennas.  The  Green’s  functions  found  in  this  paper  along  with 
the  MPIEs  can  be  applied  to  the  study  of  any  two  dielectric  layer  structure.  The  Green’s 
functions  can  be  used  in  scattering  analysis  of  two  layer  structures,  to  study  two  layer  mi¬ 
crostrip  antennas  with  circular  or  even  arbitrarily  shaped  patches,  or  to  study  stacked 
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patch  antennas  with  microstrip  or  coaxial  feeds  instead  of  an  aperture  feed. 

Different  types  of  aperture  feeds  may  also  be  explored.  If  a  microstrip  antenna 
can  be  fed  through  an  aperture  by  a  microstrip,  can  it  by  fed  through  an  aperture  by  a 
rectangular  or  circular  waveguide?  What  kind  of  modes  would  be  excited  in  a 
waveguide  by  an  aperture  coupled  microstrip  antenna?  Would  a  microstrip  antenna 
(stacked  patch  or  single  patch)  fed  through  an  aperture  by  a  waveguide  have  desirable 
characteristics?  Would  it  be  practical?  As  with  any  research,  as  many  new  questions 
arise  as  are  answered. 

The  results  of  this  research  are  particularly  applicable  to  the  study  of  microstrip 
antenna  arrays.  Because  surface  wave  effects  are  included  in  the  Green’s  functions,  mu¬ 
tual  coupling  between  array  elements  can  easily  be  accounted  for  by  extending  the  mo¬ 
ments  method  analysis.  This  is  very  important,  since  microstrip  antennas  are  usually  em¬ 
ployed  as  parts  of  an  array. 

5.3  General  Observations. 

During  the  course  of  any  project,  it  is  often  observed  that  some  techniques  or 
methods  work  better  than  others.  Whether  they  work  better  simply  because  they  save 
time  or  arc  easier  to  use,  it  is  beneficial  to  pass  along  such  information  so  future  research¬ 
ers  may  profit  from  the  experiences  of  others.  Great  use  of  commercial  software  was 
made  during  this  research  and  writing  custom  software  was  always  a  last  resort.  Finding 
a  misplaced  comma  or  errant  semicolon  can  often  consume  several  hours  or  even  days,  as 
anyone  who  has  ever  written  software  can  attest.  The  IMSL  subroutines  were  found  to 
work  very  well,  both  when  directly  integrating  over  some  interval  and  when  used  as  part 
of  custom  subroutines.  The  documentation  was  well  written  with  clear  examples  of  how 
to  use  the  subroutines. 
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The  new  software  package  Mathematica  was  used  extensively  during  this  re¬ 
search.  The  construction  of  the  Green’s  functions,  staning  from  the  general  solutions  of 
the  vector  potentials  and  boundary  conditions  through  to  the  final  expressions,  was  ac¬ 
complished  totally  with  Mathematica’ s  symbolic  equation  solving  abilities.  Deriving  the 
Green’s  functions  with  pencil  and  paper  would  have  been  a  horrendous  task,  considering 
the  number  of  equations  that  had  to  be  derived  and  the  complexity  of  the  equations.  Typ¬ 
ical  mistakes  such  as  forgetting  to  change  a  sign  or  leaving  off  a  subscript  were  totally 
eliminated.  Once  a  derivation  process  had  been  defined  in  Mathematica  for  one  case  (a 
HED  on  interface  2b  for  example),  the  same  solution  process  could  be  used  for  all  similar 
cases  (like  a  HED  on  interface  3b).  Constructing  all  the  Green’s  functions  for  this  re¬ 
search  was  greatly  eased  using  this  technique.  Also,  functions  that  must  be  called  from 
external  libraries  in  Fortran  (such  as  Bessel  functions)  are  built  into  Mathematica. 
Mathematica’ s  ability  to  plot  functions  was  very  helpful  in  determining  the  general 
behavior  of  the  integrands  of  the  Green’s  functions.  Data  generated  by  Fortran  programs 
could  be  read  by  Mathematica  and  plotted  or  used  as  data  in  further  calculations  within 
Mathematica. 

It  was  hoped  at  the  beginning  of  this  research  to  use  Mathematica  for  all  numeri¬ 
cal  problem  solving  and  avoid  using  Fortran  at  all,  but  this  was  not  possible.  While 
Mathematica’ s  numerical  integration  function  is  very  sophisticated  (no  need  to  separate 
functions  into  real  and  imaginary  parts,  integration  paths  can  be  defined  in  the  complex 
plane,  etc.)  it  is  rather  slow.  Integrations  that  took  seconds  to  do  in  Fortran  on  a  main¬ 
frame  took  several  minutes  with  Mathematica  running  on  a  personal  computer  or 
workstation.  While  it  was  certainly  possible  to  numerically  evaluate  all  the  integrals  in 
the  Green’s  functions  for  all  the  required  values  of  R  with  Mathematica,  it  would  have 
taken  days.  Mathematica  was  very  useful  in  prototyping  many  of  the  subroutines  even¬ 
tually  written  into  Fortran.  Various  numerical  problem  solving  concepts  could  be  quick- 
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ly  explored  in  Mathematica,  rather  than  go  through  the  normal  write,  compile,  run, 
debug  cycle  of  writing  Fortran  code.  While  Mathematica  may  not  be  able  to  replace 
Fortran,  it  is  certainly  a  useful  problem  solving  tool  to  use  in  conjunction  with  Fortran. 

The  approach  to  analyzing  the  aperture  fed  stacked-patch  microstrip  antenna  in 
this  paper  was  based  on  the  MPIE  and  the  Green’s  functions  solved  in  the  spatial  domain. 
Obviously,  this  is  not  the  only  approach.  There  are  several  different  approaches,  some  of 
which  are  discussed  in  the  literature  review  in  chapter  II.  Just  as  the  MTIE  solution  gives 
certain  insights  to  this  problem  and  has  certain  advantages,  these  other  methods  may 
provide  different,  but  advantageous  points  of  view.  These  different  advantages  should  be 
kept  in  mind  as  this  line  of  research  continues. 
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The  techniques  needed  to  derive  the  Green’s  functions  for  a  complicated 
structure,  such  as  an  aperture  fed  stacked-patch  microstrip  antenna,  are  more  easily  un¬ 
derstood  if  a  simple  example  is  studied  first.  The  Green’s  function  for  a  single  HED  in  a 
homogeneous  medium,  for  which  an  analytical  solution  in  spherical  coordinates  is  avail¬ 
able,  will  be  developed  first.  The  result  will  then  be  adapted  to  cylindrical  coordinates 
which  are  better  suited  to  the  study  of  microstrip  structures. 

A.l  HED  in  a  Homogeneous  Medium 

A  unit  electric  dipole  has  an  infinitesimal  length  dx  in  which  circulates  an  electric 
current  /.  These  two  quantities  produce  a  dipole  moment  Idx  which  is  set  to  unity  (1  A 
m).  The  dipole  is  set  at  the  origin  of  a  rectangular  coordinate  system  with  the  current 
flowing  in  the  positive  j:-direction  (see  Figure  A  -  1). 
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The  volume  current  density  for  this  dipole  is 

J  =  x6(r)  (A-1) 

where  5  is  the  Dirac  delta  function  and  the  unit  dipole  moment  {Idx  =  1  A  m)  is  implied 
in  the  light-hand  side  of  (A  -  1)  for  dimensional  consistency  [12].  The  surface  current  in 
the  z  =  0  plane  for  a  microstrip  structure  is  defined  as 


I 


Js=  jdz 


(A -2) 


Using  the  cylindrical  coordinate  representation  of  the  Dirac  function 


results  in 


2;tp 


_^5(p) 


(A -3) 


Js  =x 


2np  (A  -  4) 

In  a  homogeneous  medium,  the  magnetic  vector  potential  created  by  an  electric 
current  is  parallel  to  the  current  [12]  so  that 

A  =  X  Ax  (A  -  5) 

As  will  be  seen  in  Appendix  B,  the  A^^  component  must  be  a  solution  of  the  Helmholtz 
equation 

(v^  +  k^)  Ax  =  0  (A  -  6) 

It  must  also  satisfy  the  boundary  condition 


1  /3Ax^.  dAx.\_  5(p) 

p  \  az  dz  I  2np  (A  -  7) 

where  A^^^^  is  just  above  the  z  =  0  plane  and  is  just  below  the  plane.  The  A^^  compo¬ 
nent  must  also  satisfy  the  Sommerfeld  radiation  condition. 

From  symmetry  considerations,  A^^  can  only  be  a  function  of  the  distance  r  from 
the  origin  [12].  The  analytical  integration  of  the  Helmholtz  equation  produces  the  result 
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A  -  expi-jig-) 

*471  r  (A-8) 

Using  the  Lx»rentz  gauge  condition,  the  scalar  potential  is 

V  =  —4 —  p  cos  (|>  ^  exp(-jkr) 

^  ^  ’  (A -9) 

which  is  actually  the  potential  created  by  two  point  charges.  Using  the  electrostatic  rela¬ 
tionship  linking  the  point  charge  and  dipole  potentials  (3  -  93)  gives  the  potential  of  the 
point  charge 


V  =_jl_£^pUe)  = 

47cjcoe  r  jcopE  (A  -  10) 

The  potential  for  the  point  charge  could  just  as  easily  been  determined  directly  from  the 
Helmholtz  equation  without  first  evaluating  A,  but  only  for  the  homogeneous  case.  For 
an  inhomogeneous  medium,  A  must  be  found  first  as  demonstrated  in  this  simple 
example  [12]. 


A  2  Cylindrical  Coordinate  System. 

The  magnetic  vector  potential  for  the  infinite  homogeneous  case  displays  a  spher¬ 
ical  symmetry,  as  can  be  seen  by  the  vector  potential’s  sole  dependence  on  the  coordinate 
r  in  (A  -  8).  However,  the  microstrip  structure  of  this  analysis  has  circular  cylindrical 
symmetry,  where  all  the  boundary  conditions  appear  in  the  z  =  constant  planes.  The  so¬ 
lutions  obtained  here  must  then  be  expressed  in  circular  cylindrical  coordinates  to  be  use¬ 
ful.  The  solution  of  the  Helmholtz  equation  in  these  coordinates  is  [12] 

V  =  Bn{kpp)  hi{k^z)  h2(n(l>)  (A  - 1 1 ) 

where  h  j  and  h2  are  trigonometric  functions  (cos  or  sin)  or  a  linear  combination  thereof; 
n  is  an  integer  since  periodicity  in  <j>  is  required  (n=0  corresponds  to  a  <(>-independent  so¬ 
lution);  is  a  solution  of  the  Bessel  equation  (Jn.  Yn,  Hn  \  or  a  linear  combina¬ 
tion);  and  k  and  k_  are  complex  quantities  called  spectral  variables  [12]  (which  are  the 
P  ^ 

same  as  eigenvalues)  and  must  satisfy  the  condition 
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kl  +  kl  =  k^  (A- 12) 

The  present  symmetrical  case  is  independent  of  <j),  so  n-O  and 

V  =  exp(-jk2|  z|)  (A  -  1 3) 

The  Sommerfeld  radiation  condition  of  (3  -  29)  is  met  in  the  radial  direction  by  the  Han- 
kel  function  of  the  second  kind.  The  radiation  condition  is  satisfied  in  the  z  -direction  if 

lw(kj}  <  0,  Re(kz)  >0  (A  -  14) 

The  superposition  principle  applies  since  the  Helmholtz  equations  are  linear,  therefore, 
any  linear  combination  of  elementary  solutions  \|/  is  also  a  solution  [12]. 

For  an  antenna  problem,  the  domain  is  considered  infinite  and  the  spectral  vari¬ 
ables  take  on  a  continuous  variation  [20]  rather  than  the  discrete  values  of  a  bounded 
problem.  Therefore,  the  general  solution  is  an  integration  of  (A  -  13)  over  either  spectral 

variable  k  or  k_  [12]. 

P  ^ 

Integration  over  k  produces  a  Fourier  transform  where  integration  over  k 
z  P 

produces  a  Hankel  transform.  Integration  over  k^  is  better  suited  to  geometries  with 
axial  symmetry  [12],  so  the  solution  for  A^^  will  be 


A.= 


I, 


'I{kp)\j^kp,kz)  dkp 


(A  -  15) 


where  T  is  an  arbitrary  function  that  will  be  determined  later.  The  integration  path  C  is, 
in  principle,  arbitrary  between  -oo  and  «»;  but  the  nature  of  the  integrand  and  conditions 
(A  - 14)  place  certain  constraints  on  C. 

The  integrand  of  (A  - 15)  contains  the  complex  function  k^  defined  as 

kz  =  V  k^  -  kj  (A -16) 

The  conditions  (A  -  14)  arc  used  to  determine  which  branch  (Riemann  sheet)  of  the  func¬ 
tion  (A  -  16)  is  selected.  The  requirement  of  a  negative  imaginary  part  determines  a  cut¬ 
off  in  the  complex  plane  for  k^.  Because  the  lower  part  of  the  k^  plane  maps  over  the 
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2 

whole  kz  plane,  the  branch  cut  is  the  positive  real  axis  (see  Figure  A  -  2b)  [12]. 

The  equation  of  the  branch  cut  is 

Im(k2)  =  0.  Re(k^)>0  (A  -  17) 

Using  the  notation 

kp  =  ^+jv  (A -18) 

the  corresponding  branch  cuts  in  the  k^  plane  are  found  from  (A  -  12)  to  be 

X  =  0  Vv 

V  =  0  for|  A.|  <  k  (A  -  19) 

They  have  their  beginning  at  the  ramification  points  k^  =  ±k  as  shown  in  Figure  A  -  2c. 
The  second  condition  on  k^  demands  the  real  part  to  be  positive,  leading  to  Imlkj)  <  0 
and  then  to  X,\)  >  0  [12]. 

The  forbidden  regions  in  the  three  complex  planes  are  shown  by  cross-hatching  in 
Figure  A  -  2.  These  forbidden  regions  require  the  contour  C  to  cross  quadrants  1  and  3, 
passing  through  the  origin  and  closing  upon  itself  at  infinity  (Figure  A  -  2c).  If  no  singu¬ 
larities  other  than  the  ramification  points  ±k  arc  present,  the  contour  can  be  deformed 
into  the  real  axis  X,  [12]. 


■0 


(a)  (b)  (c) 

Figure  A  •  2  Transformation  of  complex  planes  with  existing  branch  cuts  and  forbidden 

2 

regions;  (a)  plane,  (b)  kz  plane,  (c)  k^  plane  [i2|. 
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It  must  also  be  noted  that  the  Hankel  function  creates  an  additional  branch  cut  in 
the  plane  located  on  the  negative  real  axis.  This  can  easily  be  removed  as  will  be 
shown  later. 

The  function  T  can  be  found  now  that  the  topology  of  the  k^  plane  has  been  de¬ 
termined.  Applying  the  uniqueness  theorem  by  identifying  (A  - 15)  with  the  previously 
obtained  solution  (A  -  8)  yields 


|i  exp(-jkr) 
47t  ^ 


I 


T^kp)  Ho  \kpp)  exp(-jk2l  z|)  dkp 


(A  -  20) 


The  value  of  T  can  be  found  by  using  a  Hankel  transformation  [12],  producing  for  the 


vector  potential 


Ax 


^  z!)  = 

j*^z 


p  exp(-jkr) 
47t  r 


(A -21) 


The  integrand  of  (A  -  21)  becomes  singular  at  the  points  k^  =  ±k,  but  these  singu¬ 
larities  are  integrable  (zero  residue).  The  singularity  of  the  Hankel  function  is  eliminated 
by  the  factor  k^  in  the  integrand,  therefore,  the  contour  C  can  be  placed  along  the  real 
axis  k  =>,[12].  Noting  this  and  using  the  relationship 


I  H?^Xp)f(?c^)x"^'dX  =  2|  J„(Xp)f( 


,  2l  n+1  , 

>.  X  dX 


produces  a  form  better  adapted  to  numerical  calculations. 


(A -22) 


A.  =  Jo(Xp)  i  exrt-ul  x|)  dX  =  i 

Jo  (A  -  23) 

where 

u=jkz=  VkJ-k^  (A -24) 

This  last  expression  for  A^^  greatly  eases  numerical  integration  procedures  by  re¬ 
ducing  the  integration  to  a  semi-infinite  interval.  The  physical  significance  of  these 
developments  is  that  a  spherical  wave  can  be  expressed  as  an  infinite  sum  (integral)  of 
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cylindrical  waves  which  have  transverse  wave  numbers  (eigenvalues)  k  varying  contin- 

P 

uously  from  0  to  <»  [12]. 


A. 3  H ED  in  a  Microstrip  Structure. 

There  are  two  main  differences  between  the  homogeneous  case  discussed  above 
and  the  case  of  a  HED  in  a  stratified  medium.  The  first  is  the  vector  potential  is  no  long¬ 
er  parallel  to  the  dipole  [20].  A  vector  potential  only  parallel  to  a  boundary  is  insufficient 
to  simultaneously  meet  the  boundary  conditions  of  continuous  A|j  and  V  across  the 
boundary  (A||  is  the  component  of  A  parallel  to  the  boundary).  These  boundary  condi¬ 
tions  can  only  be  met  if  a  second  component  of  the  vector  potential  perpendicular  to  the 
boundary  is  present.  The  second  difference  is  the  Sommerfeld  radiation  condition  no 
longer  applies  for  the  z-dependence  within  the  finite  thickness  dielectric  layers.  In  the 
infinite  medium  the  dependence  is  still  of  the  form  exp(-uz),  but  within  the  dielectric 
layers  the  dependencies  are  of  the  form  sinh(uz),  cosh(uz)  or  a  linear  combination. 
Which  function  (sinh,  cosh  or  both)  to  use  is  dependent  on  the  boundary  conditions  (see 
Appendix  B). 
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For  a  time  dependence  of  Maxwell’s  equations  take  the  form  [10] 


VxE  =  -ja)pH  (B-i) 

VxH=j03eE  (B.2) 

V  -  E  =  0  (B-3) 

V  -  H=0  (B-4) 


for  the  source  free  regions  of  the  microstrip  antenna.  The  sources  only  exist  at  the 
boundaries,  therefore,  their  effects  only  appear  in  the  boundary  conditions.  Since  both 
electric  and  magnetic  sources  must  be  accounted  for,  they  will  be  discussed  separately. 
The  boundary  conditions  for  electric  sources  will  be  developed  first  followed  by  the 
boundary  conditions  for  magnetic  sources. 

B.l  Boundary  Conditions  When  Only  Electric  Sources  Are  Present. 

The  tangential  E  and  H  fields  at  a  boundary  must  be  continuous  unless  a  source 
(electric  or  magnetic)  is  present  [10].  If  only  an  electric  source  is  present,  the  tangential 
magnetic  field  must  be  discontinuous  by  the  value  of  the  surface  current  on  the  boundary. 
For  the  case  of  a  surface  current  at  interface  2b  and  no  surface  current  at  interface  3b,  the 
boundary  conditions  are 

zx(E^2-Ey  =  0  (B-S) 

zx(H52-Hy  =  J2  (B-6) 

evaluated  at  z  =  b^j^  and 

zx(E52-Ey  =  0  (B-7) 

zx(H52-Hy  =  0  (B-g) 

evaluated  at  ^  =  b2jj. 
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The  tangential  magnetic  field  is  continuous  at  interface  3b  w.iere  there  is  no  surface 
current  present.  Since  there  is  no  electric  field  in  a  prefect  conductor,  the  tangential 
electric  field  at  the  ground  plane  (interface  lb)  must  be  zero; 

2  X  Ei2  =  0  (B  -  9) 

An  equation  similar  to  (B  -  6)  could  be  used  to  determine  the  electric  current  on 
the  ground  plane,  but  (B  -  9)  completely  defines  the  perfect  electric  conductor  boundary 
[12].  There  is  no  need  to  find  the  electric  surface  current  induced  on  the  ground  plane  for 
this  analysis. 

Maxwell’s  equations  (B  -  1)  -  (B  -  4)  along  with  the  boundary  conditions  for  the 
tangential  fields  (B  -  5)  -  (B  -  9)  and  the  continuity  equation  (3  -  13)  completely  define 
the  problem  of  the  diffracted  fields.  The  boundary  conditions  for  the  normal  fields  arc 
derived  from  the  above  equations  by  taking  the  transverse  divergence  of  (B  -  5)  -  (B  -  9) 
[12]  to  obtain 


2  •  (e2b  E22  ••  Elb  E12)  =  <l2 

(B  -  10) 

Z  (li2bH52-mbH?2)  =  0 

(B-  11) 

Z  (e3bE^2-£2bEy  =  0 

(B  -  12) 

2  ■  (p3b  H32  -  }l2l,  H22)  =  0 

(B-  13) 

at  interfaces  2b  and  3b,  and 

2-  H^2  =  0 

(B-14) 

at  the  ground  plane. 

The  corresponding  boundary  equations  for  the  magnetic  vector  potential 

arc  now 

found  by  using  (3  -  5),  (3  -  7)  and  the  Lorcntz  gauge  condition  (3  -  91) 

E  =  -jouA  -  VV 
H=-LVx  A 

V-  A  +  joieiiV  =0 

I 
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which  are  repeated  here  for  the  convenience  of  the  reader.  Using  these  equations  with 
Maxwell’s  equations,  (B  - 1)  -  (B  -  4),  it  is  seen  that  the  vector  and  scalar  potentials  are 
solutions  of  two  homogeneous  Helmholtz  equations  [12]: 

(v^  +  k^)A  =  0  (B-15) 

(v^  +  k^)v  =  0  (B-16) 

where  k  is  the  wave  number  of  the  respective  medium. 

By  applying  (3  -  5)  to  (B  -  9) 

z  X  (jo)Ai2  +  VV^2)  =  0  (B  -  17) 

is  obtained  at  the  ground  plane.  Since  U  in  (3  -  5)  is  defined  within  in  arbitrary  constant, 
the  value  V^2  =  0  is  assigned  at  the  ground  plane  and  (B  - 17)  becomes 

V^2  =  0  (B-18) 

z  X  Ai2  =  0  (B  -  19) 

at  z  =  0.  Equation  (B  - 19)  limits  (3-15)  to  having  only  the  sinh(ujjj  z)  term,  since  the 
cosh(u2jj  z)  term  does  not  go  to  zero  at  the  ground  plane  and  violates  (B  -  19).  Using 
(B-  18)  with  the  Lorcntz  gauge  condition  determines  the  boundary  condition  for  the 
normal  component  of  A12  at  the  ground  plane: 

_  p 

(B  -  20) 

which  implies  the  cosh(u2jj  z)  term  in  (3  - 16). 

At  interface  2b  (z  =  bjjj),  (3  -  5)  is  used  in  (B  -  5)  to  obtain 

zx(-jcoA22-VV^+j©Ai2  + V\^2)  =  0  (B-21) 

which  expands  to 

jtoA^2  +  —  2  +  -^  V^2 

dx  dx 

ja)A^2  +  —  =  jo>Ayi2  +  -^  V12 

5y  (B  -  22) 
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The  work  required  to  move  a  point  charge  through  the  conservative  part  of  the  electric 
field  from  infinity  to  either  side  of  the  boundary  must  be  the  same  [21].  Therefore,  the 
scalar  potential  must  be  continuous  across  the  boundary  and  at  interface  2b  (B  -  22) 
becomes 

Similarly  at  interface  3b 


Applying  (3  -  7)  to  (B  -  6)  yields 

( _1_  _  _1_  ^Ax12^  _  /  1  3A^2  _  ^Az12\  _  _  j 

\|i2b  |2ib  /  lM-2b  pib  dx  /  (B  -  27) 

and 

/  1  3Ay22  1  3Ayi2|  /  1  8A^2  1  ^Azi2j_Q 

l^i2b  4ib  ^z  l\^2b  '^ib  dx  I  (B-28) 

where  the  left-hand-side  of  (B  -  27)  is  the  y  -directed  magnetic  field  due  to  an  x  -directed 
electric  current.  With  the  current  in  the  negative  x-direction,  the  magnetic  field  points  in 
the  positive  y  -direction,  obeying  the  right-hand  rule  convention.  The  left-hand-side  of 
(B  -  28)  is  the  x  -directed  magnetic  field  and  is  equal  to  zero  since  the  coordinate  system 
is  oriented  along  the  x  -directed  current.  Since  only  two  components  of  the  magnetic 
vector  potential  are  needed  to  completely  define  the  electromagnetic  fields  [20],  Ay  is  as¬ 
sumed  uniformly  zero  leaving 

[  1  3A^2  1  3Az12\_q 

\H2b  dx  nib  dx  I  (B-29) 

which  reduces  (B  -  27)  to 

1  ^A^2  1  ^A^i2  _  j 

lA2b  dz  ^lb  dz  (B  -  30) 


V22  =  Vl2 

(B  -  23) 

.  btan  .  btan 

A22  =  Ai2 

(B  -  24) 

V^2  =  V^ 

(B  -  25) 

.  btan  .  btan 

A32  =  A22 

(B  -  26) 
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at  interface  2b.  Using  (3  -  7)  on  (B  -  8)  similarly  yields 

1  _  1  ^Ax22 

|i3b  3z  ~jX2b  5z  (B-31) 

at  interface  3b.  Finally,  applying  the  Lorentz  gauge  condition  to  (B  -  23)  and  (B  -  25) 
produces 

—L^  V-  A?2  =  V-  A^2 

eib  Vi\b  e2b  H2b  (B  -  32) 

from  which  (3  -  24)  is  derived,  and 

— A22  =  — ^ —  V-  A22 

e2b  \i2b  £2b  li2b  (B  -  33) 

from  which  (3  -  28)  is  derived. 

The  scalar  and  vector  potentials  are  related  to  the  surface  charge  distribution  at 
interface  2b  by  using  (3  -  5)  with  (B  - 10)  to  get 

jCD(e2b  A^2  *  £lb  Az12)  + e2b Elb -q2 

oz  oz  (B  -  34) 

The  boundary  conditions  for  electric  sources  at  the  other  interfaces  are  found  similarly. 

The  two  Helmholtz  equations  (B  - 15)  and  (B  -  16)  along  with  the  above  derived 
boundary  conditions  define,  in  principle,  both  potentials.  However,  in  order  to  obtain  a 
unique  solution  an  additional  constraint  must  be  imposed.  This  constraint  is  known  as 
the  Sommerfeld  radiation  condition 

where  y  is  any  scalar  solution  of  Helmholtz’s  equation.  The  radiation  condition 
constrains  the  solutions  to  those  that  propagate  away  from  the  sources  and  decrease  with 
distance  [12], 

From  these  boundary  equations,  it  may  appear  the  scalar  potential  V  is  associated 
with  the  surface  charge  distribution  q  and  the  vector  potential  A  is  associated  with  the 


jky  =0 
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surface  current  J.  But,  where  the  vector  potential  A  is  determined  solely  in  terms  of  the 
surface  current  J,  the  scalar  potential  V  cannot  be  expressed  only  in  terms  of  the  the 
surface  charge.  Because  of  the  presence  of  the  term  in  (B  -  34)  the  determination  of  V 
requires  previous  knowledge  of  the  vector  potential  A  normal  to  the  currents.  In  a 
homogeneous  medium  where  the  permittivity  and  permeability  are  not  discontinuous,  the 
normal  component  A^  vanishes  and  V  is  independent  of  A.  In  an  inhomogeneous 
structure,  the  normal  component  A^  is  always  present  [12]. 

B2  Boundary  Conditions  When  Only  Magnetic  Sources  Are  Present. 

For  the  case  when  only  a  magnetic  source  is  present,  the  tangential  electric  field 
must  be  discontinuous  by  the  value  of  the  magnetic  surface  current  on  the  boundary.  The 
only  magnetic  sources  for  this  analysis  are  the  equivalent  magnetic  sources  used  to  model 
the  aperture  in  the  ground  plane.  Since  the  magnetic  sources  are  on  the  ground  plane, 
special  care  must  be  taken  in  deriving  the  boundary  conditions  at  interfaces  lb  and  la. 
The  boundary  condition  at  interface  lb  is 

-  z  X  Ei2  =  Ml  (B  -  35) 

The  corresponding  boundary  equations  for  the  electric  vector  potential  are  found 
by  using  (3  -  6),  (3  -  8)  and  the  Lorentz  gauge  condition  for  magnetic  sources  (3  -  92) 

E=-lVxF 

e 

H  =  -jtoF  -  VVn, 

V-  F  +jo)e^iVm  =  0 

repeated  here  for  convenience.  Separating  (B  -  35)  into  x  and  y  components  produces 

e5i2  =  (B  -  36) 

E$i2  =  0  (B-37) 
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where  the  left-hand  rule  convention  for  magnetic  sources  is  demonstrated  in  (B  -  36). 
Applying  (3  -  6)  to  these  two  equations  yields 


^aF$i2 


eib  \  3z  d\  } 


xl 


(B  -  38) 


and 


1  (3Fz12  ^Fyi2\_Q 
eib  \  3x  3z  / 


(B  -  39) 

Only  two  components  of  the  electric  vector  potential  are  needed  to  completely 
define  the  electromagnetic  fields,  therefore,  Fy  is  assumed  uniformly  zero  leaving 


zl 

3x 


*ai=o 


(B-40) 


which  reduces  (B  -  38)  to 


eib  oz 


(B.41) 

at  z  =  0.  The  z-dependence  of  F  in  the  dielectric  layers  for  this  problem  is  in  a 
cosh(ujjjZ)  or  sinh(ujjjZ)  term  which  is  not  affected  by  the  d/dx  operation  in  (B  -  40), 
therefore,  only  the  sinh(u2jjZ)  term  can  be  used  in  (3  -  61)  for  F212  to  meet  the  boundary 
condition.  Also  notice,  since  the  source  is  on  the  ground  plane,  the  magnetic  scalar 
potential  cannot  be  set  to  zero  at  the  ground  plane  as  with  the  electric  sources  at  the 
dielectric  boundaries.  The  boundary  conditions  at  interface  la  are  the  same  as  above  and 
the  boundary  conditions  for  the  dielectric  boundaries  (2a,  2b,  and  3b)  are  found  from  the 
boundary  conditions  for  electric  sources  by  using  duality. 
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C.l  HMD  on  Interface  la. 

N3xi(>.)  =  Eia  X  exp(bia  U2a)  Sech(bia  Uia) 

N^,(X)  =  eia  X  [Ula  +  eal2  U2a  tanh(bia  Ula)] 

N^iW  =  2eia>.^l  -|^al2eal2)cSCh(2biaUia) 

N^i(x)  =  e2a  exp(biaU2a)(l  *  ^al2  Can)  Sech(bia  Uu) 

C.l  HED  on  Interface  2a. 

Ntx2W  =  ^lla  X  csch(bia  Uu) 

Ndx2(^)  =  ^^la  exp(bia  U2a) 

NL2(>^)  =  JXla  (1  -  |Ial2  eal2)  Sech(bia  Uia) 

Ndz2(^)  =  H2a  X^  exp(bia  U2a)(l  -  Hal2  Ea^) 

C.3  HMD  on  Interface  lb. 

Nb<lW  =  -  Eib  X  [eb23  U3b  cosh(b2b  U2b)  +  U2b  sinh(b2b  U2b)]  sech(bib  uib) 
nJ„(x)  =  Eib  X  [U2b  cosh(b2b  U2b)  +  eb23  U3b  sinh(b2b  U2b)]  sech(bib  uib) 
nLiW  =  Eib  X  U2b  exp(b2b  U3b)  sech(bib  uib) 

j^b  =  ej^  X.  (  tanh(bib  uib)]  U2b  cosh{u2b  (b2b  -  bib))  1 

1+ [eb23  uib  U3b  +  ebi2  u^b  tanh(bib  uib)]  sinh{u2b  (b2b  -  bib))| 

|[-  M'bn  ebi3  +  cosh^{u2b  (b2b  -  bib))]  2  uib  csch(2  bib  uib) 

2  I  +  [-  M-bl2  £bl2  uib  -  ^bl3  ebl3  U^b  +  Hb23  eb23  U^b] 

nLi(>.)  =  eib  X,  1*2  csch(2  bib  uib)  sinh^(u2b  (b2b  *  bib)) 

+  [M-b23  -  ^bl3  ebl2  -  M-bl2  EbO  +  eb23] 

*  U2b  U3b  csch(2  bib  Uib)  smh(2  U2b  (b2b  -  bib)) 
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0  -  I^b23  eb23]  Elb  Uib  U2b  csch(bib  Uib)  cosh{bib  U2b)  ' 
+  [-  I^bi2  +  ^bi3  eb23]  Eib  uib  sech(bib  Uib)  sinh{bib  U2b) 

^  (4b  13  Elb  -  4b23  e2b)  U3b  COSh(b2b  U2b) 

+  (4bl2  £lb  -  E2b)  U2b  Sinh(b2b  U2b)  , 

*  U2b  sech(bib  Uib)  cosh(u2b  {b2b  -  bib)) 

^  (4bl3  Ebn  -  4b23  eb23)  U3b  COSh(b2b  U2b) 

+  (4bl2  ebl3  -  eb23)  U2b  sinh(b2b  U2b) 
e2b  U3b  sech(bib  uib)  sinh(u2b  (b2b  -  bib))  j 


[-  1  +  4b23  eb23]  Eib  Uib  U2b  csch{bib  Uib)  sinh{bib  U2b)' 

+  [4bi2  -  4bi3  €b23]eib  uib  sech(bib  uib)cosh(bib  U2b) 

Nc  1(5^)  =  +  ^'  U2b) 

“  /  L  +  (-  4bl2  Elb  +  e2b)  U2b  COSh(b2b  U2b)  . 

*  U2b  sech(bib  uib)  cosh(u2b  (b2b  -  bib)) 

^  {-  4bl3  ebl3  +  4b23  eb23)u3b  sinh(b2b  U2b) 

+  (- 4bl2  ebl3  +  eb23)  U2b  COSh(b2b  U2b) 

\*  e2b  U3b  sech(bib  uib)  sinl^U2b  (b2b  -  bib)) 


[-  4bi3  Eib  +  e3b]  U2b  sech(bib  uib)  cosh(u2b  {b2b  -  bib))' 
^  [(-  4b23  eib  +  e3b  ebi2)  uib  csch(bib  uib) 

+  (-  4bi2  Eib  +  e2b)  U3b  sech(bib  uib)  . 

*sinh(u2b  (b2b  -  bib))  , 


C.4  HED  on  Interface  2b. 

Nm2(^)  =  4ib  ^  [  U2b  cosh(u2b  {b2b  -  bib))  +  4b23  U3b  sinh(u2b  (b2b  -  bib))]  csch{bib  uib) 
^1x2^  =  -  |Ilb  X  [4b23  U3b  C0Sh(b2b  U2b)  +  U2b  sinh(b2b  U2b)] 

=  4lb  ^  [U2b  C0Sh(b2b  U2b)  +  4b23  U3b  sinll(b2b  U2b)] 
nL(x)  =  ^ibXexp(b2b  U3b)u2b 

[1  -  4bi3  ebi3]  uib  sech(bib  uib) 

h)  M  ■‘■[(1  -  4bl2  Ebu)  uib +{- 4bl3  EbO  +  4b23  eb23)ulb] 

*  sech(bib  Uib)  sinh  (u2b(b2b  *  bib)) 

+  [4b23  -  4bl3  Ebl2 ' 4bl2  Ebl3  +  Eb23] 

*  H21di2l2.  sech(bib  uib)  sinh(2  U2b(b2b  -  bib)) 
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\ 


[-  £bl2  +  4b23  Ebis]  l^lb  U^b  sinh(bib  U2b) 

'1  -  |j,b23  eb23]  M-Ib  Uib  U2b  cosh(bib  U2b)tanh(bib  uib) 

(M-lb  ebl3  -  M-2b  £b23)  U3b  COSh(b2b  U2b) 

.  +  (-  ^2b  +  ^Alb  £bl2)  U2b  Sinh{b2b  U2b) . 

*  U2b  COSh(U2b(b2b  -  bib)) 

(^bl3  £bl3  -  Hb23  eb23)U3b  COSh(b2b  U2b) 

.  +  (-  M-b23  +  M-bl3  Ebnj  U2b  Sinh(b2b  'J2b) . 

*  ^2b  U3b  sinh(u2b(b2b  -  bib)) 


[Ebl2  -  M-b23  Ebn]  4lb  U^b  COSh(bib  U2b) 

1  +  ^b23  eb23]  V^xh  Uib  U2b  sinh(bib  U2b)  tanh(bib  uib) 

(-  ^ilb  £bl3  +  42b  £b23)  U3b  Sinh(b2b  U2b) 

+  (42b  -  4lb  £bl2)  U2b  COSh{b2b  U2b)  . 

*  U2b  cosh(u2b(b2b  -  bib)) 

^  {-  4bl3  ebl3  +  4b23  eb23)  U3b  sinh(b2b  U2b) 

+  {4b23  -  4bl3  Ebn)  U2b  C0Sh(b2b  U2b)  . 

*42b  U3b  sinh(u2b(b2b  -  bib)) 


^  |[43b  -  4lb  Ebia]  U2b  COSh(u2b  (b2b  -  bib)) 
=  exp(b2b  U3b)  U2b  ^  |+ 


(42b  -  4lbebl2)U3b  ll 

L+  (43b  4bi2  -  4ib  eb23)  Uib  tanh(bib  uib)J ( 
*sinh(u2b  (b2b  -  bib))  I 


C.J  H ED  on  Interface  3b. 


N«3(x)  =  4ib^U2b  csch(bibUib) 


nSx3(>.)  =  42b  ^[uibCosh{bib  U2b)coth(bib  Uib)-  4bi2  U2b  sinh(bib  U2b)] 


Nc^3(^)  =  42b  ^  [-  Uib  sinh(bib  U2b)  coth(bib  uib)  +  4bi2  U2b  cosh(bib  U2b)] 

nU>^)  =  |i,2b  ^  exp(b2b  U3b) 

Mb  ^ )  ..  1 2  P  ■  ^bi3  Ebn]  U2b  cosh{u2b(b2b  -  bib))  secKbib  uib) 

Naz3lA,;-^ibU2b^  /  r(ebl2-4b23  ebl3)uibCSCh(bibUib)  1 . 


4b  12  U2b  COSh(u2b(b2b  -  bib)) 

+  Uib  coth(bib  Uib)  sinh(u2b(b2b  -  bib)). 


L+  (-  4bi2  ebi3  +  eb23)  U3b  sech(bib  uib). 


1 


sinh(u2b(b2b  -  bit,})j 


-  no- 


Appendix  C:  Vector  Potential  Parameters 


/[{l^lb  EbO  -  M-2b  eb23)  M-3b  COsh(b2b  U2b)] 

L+  (-  l^2b  +  ^Ib  Ebn)  |^2b  sinh(b2b  U2b)  J 

2  +  (■  £bl2  +  lJ-b23  Ebis)  U2b  sinh(bib  U2b) 

Nbz3W  =  X  L+(l -|^b23eb23)uibCosh(bibU2b)tanh(bibUib). 

*  4ib  U2b  cosh{u2b(b2b  -  bib)) 

^  (1  -  )^b23  eb23)  Ulb  COSh(bib  U2b) 

[+  (-  £bi2  +  M-b23  Ebn)  U2b  coth(bib  Uib)  sinh(bib  U2b). 

\*  ^i2b  Ulb  Sinh(u2b(b2b  -  bib)) 

'  (-  M-lb  ebl3  +  ^2b  eb23)  ^3b  sinh{b2b  U2b)l 
L+(42b-4lbebl2)42bCOSh(b2bU2b)  J 
2  +  (Ebl2  -  |ib23  Ebn)  U2b  COSh(bib  U2b) 

NczsW  =  ^  /  L+  {- 1  +  Ub23  eb23)  Ulb  sinh(bib  U2b)  tanh(bib  uib) 

*  l^lb  U2b  COSh(U2b(b2b  -  bib)) 

+  [(- 1  +  Ub23  eb23)  Ulb  sinh{bib  U2b) 

L+  (ebl2  -  I^b23  Ebn)  U2b  COth(bib  Ulb)  COSh(bib  U2b). 

\*  |X2b  Ulb  sinh(u2b(b2b  *  bib)) 

[-  Ulb  +  ulb]  |i3b 

^  +  .M-3b  Ulb  -  M-lb  ebl3  ulb]  C0Sh^(u2b  {b2b  -  bib)) 

Ndzsl^)  =  exp(b2b  U3b)  "^  !+  •  \i-2b  eb23  Ulb  +  M-3b  M-bi2  ebi2  ulb]  sinh^(u2b  (b2b  -  bib)) 

+  [(M-3b  ebl2  *  l^2b  Ebis)  COth(bib  Uib) 

L+  (43b  4bi2  *  4ib  eb23)  tanh(bib  uib). 

*  iLL^}2b  5 U2b  (pzb  -  bib)) 


-  Ill  - 
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Dix)  =  ^al2  U2a  +  Uu  COth(bia  Uu) 

=  eal2  U2a  +  Uu  tanh(bia  Ula) 

Ula  =  V>.^-ku,  U2a  =  V 

Db(x)  =  t^‘’13U3b  +  Uib  COth(bibUib)]  U2b  COSh{u2b(b2b  -  bib)) 

+  [^bi2ulb  +  Ub23UibU3b  coth(bibUib)]  Sinh{u2b(b2b  -  bu)) 

D^>.)  =  I^^bi3U3b  +  uib  tanh(bibUib)]  U2b  cosh{u2b(b2b  -  bu)) 

+  [ebi2uib  +  eb23UibU3b  tanh(bibUib)l  sinh(u2b(b2b  *  bu)) 


uib  =  V  X^  -  k?b ,  U2b  = 


-  k2b ,  U3b  =  V  -  k|b 


D.l  Green’s  Functions  for  Electric  Fields  From  Electric  Sources. 
D.l.J  HED  at  interface  2a: 


Ga22(R)  =  ^ 

dJa.) 


GaY2(R)  =  G?S2(R) 


G522(R)  =  :r^  H?\XR)  X  H_«12  ^laBnlibuJila)]  dX 

^"^2.  dWdUx) 


For  large  X,  these  integrands  decay  as 

D.l. 1.1  Asymptotic  forms  for  small  R: 

G5a(R)  =  ^j|  lX..(RA)dX-j^| 
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Gq22(R)  =  7-*— I  f  I^22(RA)ca--— 1 - f  I^dx+ 4-7-1 - / 

’  47re2a^  ^  ’  l+eai2j^  R  Rl-^ean) 

where  Ia22(RA)  and  Iq22(Ri^)  arc  the  integrands  of  the  original  functions  with  Hq^^Xr) 

replaced  by  Jo(^R). 

D.1.2  HED  at  interface  2b: 

Ga^2(R)  ~  COSh(u2b  (b2b  ~  b)b))  +  p.b23  ^3b  sinh(u2b  (b2b  ~  bib)) 

^  D^(X) 


\ 


gS^(R)  =  G^*2(R) 


r 


G^22{R)  = 


47te2b 


D^xiDSix) 


eb23  uib  U3b  COSh^(u2b  (b2b  -  bib)) 

q.  M-b23  uib  U3b  -t-  (|ibl2  U^b  +  ltbl3  eb23  U^b) 

*  uib  tanh(bib  uib) 

sinh^(u2b  (b2b  -  bib)) 

uib  +  M-b23  Eb23  uib  +  (M-bl3  +  M-bl2  eb23) 

♦uibU3btanh(bibUib) 

*  sinh(2  U2b  (b2b  -  bib)) 


'  +  ltbi3  eb23  Uib  uib  tanh(bib  uib) 
For  large  X,  these  integrands  decay  as  X'^^. 


dX 


G^2(R)  =  ^| 


G«2(R)  =  0«2(R) 


GS32(R)  =  2:i 


47te3b 


/[u3b  +  ltbi3  Uib  tanh(bib  uib)] 
Hif\XR)Xu2b  *  “2b  cosh(u2b  (b2b  -  bib)) 

+  M-bl3  ebl2  uib  +  (Hb23  ’  M-bl3  Ebu)  uib 

I  +  ^bi2  Uib  U3b  tanh(bib  Uib) 

[*  sinh(u2b  (b2b  -  bib)) 


d^IxId^x) 


dX 


For  large  X,  these  integrands  decay  as  exp[-X  (b2{j  'bijj)]- 
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D. 1.2.1  Asymptotic  forms  for  small  R: 

G^(R)=i^|[  iWR.>-)ca-  { 

2k  I  1  +  libl2  +  M-b23  +  I^bl3 


Jo(x) 

R 


dx 


J _ 1  +  M-b23 _ 

R  1  +  |Abl2  +  4b23  +  l^bBl 


G$22{R)  =  f  WRA)  dX  -  - - 

27ie2b  1  +  £bl2  +  eb23  +  ^*>13  R 


dx 


1 


1  +  eb23 


R  1  +  £bl2  +  £b23  +  £bl3 


D.1.3  HED  at  interface  3b: 


G^*3(R) 


_  M-2b 
4k 


Hn'(>-R)>-  [Mbn  U2b  COSh(U2b  {b2b  -  bib)) 

[+  Uib  coth(bib  Uib)  sinh{u2b  (b2b  -  bib)). 


dX 


oS(R)  =  G^5(R) 

f 


GS33(R) 


4:ie3b 


U3b  +  ^ibi3  uib  tanl<bib  uib) 
[*cosh^(u2b  (b2b-bib)) 

(u?b  +  ^bl2  Cbl2  uib)  U3b 


uit 


(2y  I  L'*’ ^b23  ebl2  Uib  U^b  COth(bib  Uib) 

bn  (  *  (b2b  -  bib)) 

Dnl^)  1  rnb23  Uib  +  ^bl3  Ebl2  U^b 

ebi2C0th(bibUib) 


*^sinh(2  U2b  (b2b-bib)) 


For  large  X,  these  integrands  decay  as 


G^a(R)  =  Gi^(R)  =  G^*2(R) 


dX 


Appendix  D:  Complete  Green's  Functions 


gS33(R)  =  T^lf  lyRAlca-- - Ltlm - [  JoW 

2ne3b|j^  1  +  ebi2  +  eb23  +  R 

,  1  l+£bl2  1 

R  1  +  EbU  +  eb23  +  £bl3  I 

For  large  X,,  these  integrands  decay  as  exp[-X  (b2{j  -bjjj)]. 

D. 1.3.1  Asymptotic  forms  for  small  R: 


dx 


g£5(R)=^ 

ZK 


(  I^3(R,>.)dX.-- - L±i^bi2 - 1*  ^ 

jg  1  +  I^bl2  +  lib23  +  ItboJ^  R 


dx 


^  J _ ^  +  Ubl2 _ 

R  i  +  )J-bl2  +  ltb23  +  l^bB 


G533(R)=:rS— jf  i533(R,>^)ca-- - - f  ^ 

27te3b  I  1  +  £bl2  +  £b23  +  £bl3  R 


dx 


+i- 


1  gbl2 


R  1  +  ebi2  +  eb23  +  £bl3 


D.2  Green’s  Functions  for  Magnetic  Fields  From  Magnetic  Sources. 
D.2.1  HMD  at  interface  la: 

Gr*i(R)  =  ^  [  H^Xr)  X  eal2U2at^i^buUu)  ^ 


471 


Ui, 


D*J[X.) 


Gri"i(R)  =  Gril{R) 


Gmll(R)  =  -;r^ 


47tpu 


uiaD3(x)D‘j[x) 


(^.12  +  eal2)  Uu  U2,  +  U?,  COth(bu  Ula) 
+  2{-  1  +llal2eal2)uiaCSCh(2biaUia) 

.+  Iiai2  eai2  uL  tanh(bia  Ula) 


dX 


For  large  X,  these  integrands  decay  as  X’^^. 
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D.2.1.1  Asymptotic  forms  for  small  R: 


Gmll(R) 


D.2.2  HMD  at  interface  Ib: 

r 


G^pI^(R)  =  ^ 


/[uib  +  ebi3U3btanh(bibUib)] 

^  fa/  1  *  U2b  coshiu2b  (b2b  -  bib)) 
uib  DitAM  U  [eb23  uib  U3b  +  ebi2  ulb  tanh(bib  uib)]j 
\  *sinh(u2b(b2b-bib)) 


Gjn(R)  =  GFn{R) 


-  116- 


ca 


Appendix D:  Complete  Green's  Functions 


G“ii(R)  = 


_  1  /  h^^\xr]x 


U2b 


I  |Ibl3  UlbU3b  +  UibC0th(bibUlb) 

+  (-  1  +  4bi3  ebi3  )  2  u?b  csch(2  bib  uib) 

.+  4bi3  ebi3  U3b  tanh(bib  uib) 

+  fbl3  Uib  U2b  U3b  COSh^{u2b(b2b  '  bib)) 

(4bl3  +  |^b23  £bl2  +  4bl2  eb23)  Uib  U^b  U3b 
+  (u^b  +  I^b23  eb23  U^b)  Uib  COth(bib  Uib) 

+  ((-  1  +  M-bl2  ebl2)  uib  +  {^bl3  ebl3  *  |ib23  Eb23)  ulb) 

*2  Uib  csch(2  bib  Uib) 

+  (4bi2  ebi2  uib  +  4bi3  ebi3  uib)  uib  tanh(bib  uib)  J 

*  sinh^(u2b(b2b  -  bib)) 

(M’bll  +  £bl2)  Uib  uib  +  (M-b23  £bl3  +  M-bl3  Eb23)  Uib  uib 
+  (^b23  +  eb23)  Uib  U3b  COth(bib  Uib) 

+  (-  ^ib23  +  4bl3  ebl2  +  M-bl2  ebl3  -  eb23)  2  ufb  U3b  CSCh(2  bib  Uib) 
.+  (4bl3  £bl2  +  ^bl2  Ebo)  uib  U3b  tanh(bib  Uib) 
\*^sinh(2u2b(b2b-bib)) 


I 


For  large  X,  these  integrands  decay  as 
D.2.2.1  Asymptotic  forms  for  small  R  ; 

G^^^(R)  =  ^||  I^ii(R,x)dX-|  +  ^ 


G^i 


D.3  Green’s  Functions  for  Electric  Fields  From  Magnetic  Sources. 
D.3.1  HMD  at  interface  la: 

Ge2i(R,C)  =  -  sin(2C)[l  IL2i(R)  -  Ia2l(R)] 

GS(R,C)  =  -  [-  iL2i(R)  -  cos2(c)  I^^i(R)  +  cos{2C)  I*.2i(R)] 
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Gki(R.C)  =  -  [iL>i{R)  +  sin\)  fa2i(R)  +  co42C)  I^iifR)] 

Ge2^(R>C)  =  -  sin(2^)j^-  +  I^i(R)j 

where 


l'a2l(R)=  I  Hn>^R)  ^  ^  SecKbia  Uia)  ^ 

dJxId'JIx) 


Iti(R) 


^^(1  -4al2eal2)sech(biaUla) 

R  dJx|d*Jx) 


IL2i(R) 


h”'|xr)  sect<bi,U|.)  ^ 


For  large  X,  these  integrands  decay  as  exp[-X  b^^]. 

D.3.2  HMD  at  interface  lb: 

0|5i(R.?)  =  -  sin(2C)[i  lt2i(R)  -  ft2l(R)] 

Gki(R'0  =  ‘  [-  Ib2l(R)  -  cos^^)  Ib2l(R)  +  COs{2^)  Ib2l(R)] 

^Ki(R-?)  =  ■  [lb2l{R)  +  sin^^)  Ib2i(R)  +  cos{2^)  Ib2i(R)] 

Ggi;(R,C)  =  -  sin(2C)[-  r  l'b2i(R)  +  Imi(R)] 
where 
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Ib2l(R)  = 


sech(bib  uib) 

(1  -  ^bl3  £bl3)u2b 

+  [(1  -  ^ibl2  £bl2)  U2b  +  (*  4bl3  ebl3  +  ^b23  £b23)  U^b 

*  (  *  sinh^(u2b  (b2b  -  bib)) 

+  [|^b23  -  4bl3  £bl2 ' 4bl2  £bl3  +  eb23] 

*  —2^'"  sinh{2  U2b  (b2b  -  bib)) 


dk 


js  j  hH^R)  secKbib  Uib) 

RDaxlD^X) 


(1  -  4bl3  ebl3)u2b  ' 

+  [{1  -  M’bn  Ebn)  uib  +  (-  H-bl3  EblS  +  M-b23  Eb23)  U^b] 
*  sinh^{u2b  (b2b  -  bib)) 

+  [^b23  -  I^bl3  ebl2  '  M-bl2  £bl3  +  eb23] 
*^^^^sinh(2  U2b(b2b-bib)) 


dk 


,rRi=  I  hH>^r)>^  U2b  sech(bib  Uib)  rebi3  ujb  cosh(u2b  (b2b  -  b]  b))  1 

j  D^X)  +ebi2U2bsinh(u2b(b2b-bib))J^ 

For  large  X,  these  integrands  decay  as  exp[-X 


G^^(r,C)  =  -  sin(2C)[l  I[.3i{R)  -  Ib3i(R)] 

Gki(R-C)  =  -  [-  Ib3l(R)  -  cos^c)  Ib3l(R)  +  COs(2C)  Ib3l(R)] 
Gki(R.C)  =  -  [lb3i(R)  +  Sin2(;)  rb3i(R)  +  cos{2C)  ibi(R)] 

Ge3i(R»C)  =  -  sin(2^)j^-^Ib3i{R)  +  Ib3i(R)] 

where 


Ib3l(R)  = 


H 


(2 


\xr)x~ 


d^IxId^x) 


[[1  -  Hbi3  Ebn]  uib  sech(bib  uib)  cosh{u2b  {b2b  -  bib))i 
r(ebl2  -  ^b23  Gbo)  Uib  U2b  CSCh{bib  Uib)  1 

[+  (-  ^ibi2  ebi3  +  eb23)  U2b  U3b  sech(bib  uib). 

\*  sinh(u2b  {b2b  -  bib)) 


dX 
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Ib3l(R)  = 


rd^Wd^x) 


[1  -  ^bi3  Ebn]  U2b  sech(bib  uib)  cosh(u2b  (b2b  -  bib)) 
+  [(ebi2  -  l^b23  ebi3)  uib  U2b  csch(bib  uib) 

[+  (-  I^bi2  ebi3  +  eb23)  U2b  U3b  secKbib  uib). 

\*  sinh{u2b  (b2b  -  bib)) 


U3b  sech(bib  uib)]  dX 


Ib3l(R)  = 

For  large  X,  these  integrands  decay  as  exp[-X  bjjj]. 


D.4  Green’s  Functions  for  Magnetic  Fields  From  Electric  Sources. 
D.4.1  HED  at  interface  2a: 

G?fi"2(R,C)  =  sin(2C)[^  Iai2(R)  -  lii2{R)] 

Ghi2  (r,C)  =  [iLi2(R)  -  cos  ^C)l!ii2(R)  +  cos(2C)lii2(R)] 

OhuM  =  [-  Ial2(R)  +  sin^C)  IIi2(R)  +  cos{2C)  I1i2{R)] 

GK(R.C)  =  sin(2C)[- 1  I1i2(R)  +  Ial2(R)j 
where 

Ial2(R)=  {  H^i^^Xr)  — ^al2)  secb(bia  Uu)  ^ 


I|i2(R)  = 


hHxr)  'X  (1  -  |Ial2  eal2)  Sech(biaUia) 

R  d3(x)dUx) 


dX 


iLi2(R)= 

Jc 

For  large  X,  these  integrands  decay  as  exp[-X  b|^]. 


dX 


-  120- 


AppendixD:  Complete  Green's  Functions 


D.4.2  HED  at  interface  2b: 

Ghi2(R.(;)  =  sin(2C)[i  Iti2(R)  -  Ibi2(R)] 

Ghi2(R.I;)  =  [lbi2(R)  -  cos^;)  lt,2(R)  +  cos(2;)  I^,2(R)] 
GhuIR.^)  =  [-  Ibl2(R)  +  Ibl2(R)  +  COs(2C)  I^12(R)] 
Ghi2(R.?)  =  sin(2C)[-  \  'bi2(R)  +  Ih2(R)] 


where 


Ibl2(R)  = 


Ibl2(R)  = 


[1  -  }ibi3  Ebo]  U2b  sech(bib  uib) 

Hn>-R)  ebl2)  uib  +  {-  I^bl3  ebl3  +  |ib23  Eb23)  U3b] 

D^X)  *  sech(bib  uib)  sinh^(u2b(b2b  -  bib)) 

+  [M-b23  -  M-bl3  Ebl2  '  Ubl2  Ebn  +  Eb23] 

*  H2k^  sech(bib  uib)  sinh(2  U2b{b2b  -  bib)) 


|[1  -  |ibi3  Ebn]  u^b  sech(bib  uib) 

Hf  \xr)  +  [(1  -  ^bi2  ebi2)  uib  +  (-  ^bi3  ebi3  +  ^b23  eb23)  uib] 

R  D^fx)  D^X)  *  sech(bib  uib)  sinh^{u2b{b2b  -  bib)) 

+[^b23-Hbl3ebl2-^ibl2ebl3  +  Eb23] 

*  sech(bib  uib)  sinh(2  U2b(b2b  -  bib)) 


t  j  .  h?Mx  uib  [  U2b  cosh(u2b  {b2b  '  bib))  1  cscWb  u  )  dX 
”  dIx)  ^"^3  U3b  Sinh(u2b  (b2b  -  bib))J 


For  large  X,  these  integrands  decay  as  exp[-X  bjjj]. 

D.4.3  HED  at  interface  3b: 

Gh13(R.O  =  sin(2C)[^  Ibl3{R)  -  Ibl3{R)] 

Gm3  (r,;)  =  [i{,13(R)  -  cos  Wlbi3(R)  +  cos{2C)l^i3(R)] 
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Ghi3(R>C)  =  [-  Ibi3(R)  +  sin^c)  Ibi3{R)  +  cos(2C)  Ibi3(R)] 

GH73(*^’d  =  sin(2i^)|^- 1  Ibi3{R)  +  Ibi3(R)] 
where 


Ibl3(R)  = 


f 


J  C 


U2b 


-  4bi3  Ebn]  U2b  cosh(u2b{b2b  -  bib))  sech(bib  uib)' 
(ebi2  -  I^b23  Ebn)  uib  csch(bib  uib) 

.+  {-  ^bi2  ebi3  +  eb23)  U3b  sech(bib  uib) 
sinh(u2b(b2b  -  b,b)) 


dX 


r 


Ibl3(R)  =  J 


RDax)D‘j(),) 


-  M-bn  Ebn]  U2b  cosh{u2b{b2b  -  bib))  sec 
(ebl2  -  I^b23  ebl3)  Uib  CSCh{bib  Uib) 

.+  (-  M-bi2  ebi3  +  eb23)  U3b  sech(bib  uib) 

sinh(u2b(b2b  -  bib)) 


Kbib  uib)' 


/ 


dX 


i!,13(R)=  Hn)^R)?^HlbU2bCsch(bii,Uii,)^ 

For  large  X,  these  integrands  decay  as  exp[-X  bjjj]. 
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<il  |  5j„(plp').Tij(p')ds' 

+j^[  [G5„(pTiip')-G5„(piiip-)]n,ip')ds' 

J  Sii 


[  ^lllplpj  T,j(p')ds' 

Vcii  /Slj 

+ j^f  [oS.,(priip')  -  GS„(piiip')]  n.ip')  ds' 

J  Sli 


-U'l  5izi(plp')- T(j(p')ds' 

+  j  [Gq22(pfilp')  -  Gq22(pf,lp')]  Ilf/p')  ds' 

<JI'  [  Slzilplpj  T2j(p')ds' 

^^Jc2i  Js2j 

*  I  [GS“(p2iip')  -  G$22(p«lp')]  n2,(p')  ds' 

Js2i 


r=ff  d|.f  S;23(plp'^  T3j(p-)ds' 

^Vc2i  Jsij 

+  Jo)^^  I  [^S23{p2ilp')  -  GS23(p2ilp')]  n3j(p')  ds' 

/Sjj 


f  dl-  [  ^2(plp'}T2j(p')ds' 

+  j(oa\b2  I  [^S32(p3ilp')  -  Gj32(p3ilp')]  n2j(p')  ds' 
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[  ^sipip'^  Tjjip'jds' 

■’vcsi  hii 

jwk^l  [^WpJilp')  -  GS,3(pi,lp')]  nj/p'l  ds' 

Js3i 


pal  2 _ |w 

'-ij  - 


pa21 _ pw 

'-ij  - 


=  ®f  dl-  [  ^izlplp').  Tr/p')ds' 

Jcii  /s§ 

"iff  f  5j2l(plp')-Ti/p')ds' 

^WC2i  Jsij 


'''”°ifi  ‘“'1 

'iff  f  Tij(p')ds' 

Jsij 


pb21  _  jw 

'“‘j  “ 


c"/‘=pf  dl-  [  G|3,(plp')T,j(p')ds' 

'■Vc  J«, 


fC3i  »Slj 


I 

i 
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A  description  of  the  programs  and  subroutines  developed  during  this  research  are 
presented  in  this  appendix.  All  practical  numerical  integration  routines  were  accom¬ 
plished  in  Fortran  77  on  an  ELXSI  System  6400  computer  with  EMBOS  13,  running 
under  the  UNIX  operating  system.  Copies  of  the  actual  Fortran  code  may  be  obtained  by 
contacting 

Maj  Harry  Barksdale 
AFIT/EN 

Wright-Patterson  AFB,  OH 
45433 

Routines  written  in  Mathematica  for  the  graphical  display  of  data  are  also  presented  in 
this  appendix.  The  Mathematica  routines  were  developed  on  a  Macintosh  SE/30 
personal  computer,  but  should  operate  without  modification  in  Mathematica  running  on 
any  platform. 

F.l  GREENMAIN 

GREENMAIN  is  the  main  program  used  to  numerically  solve  the  integrals  in  the 
Green’s  functions  Ga22(R),  Gq22(R),  G^3{R),  Gq23(R),  Ge2ii(r,C),  and  Ge2i(R.C).  All 
material  parameters  of  the  antenna  (thickness,  permeability,  and  permittivity  of  the 
dielectric  layers,  observation  frequency,  etc.)  are  set  in  GREENMAIN.  The  critical 
point,  Xj,,  for  the  asymptotic  expressions  of  G^2(R)  and  Gq22{R)  are  found  and  the 
integrations  in  the  asymptotic  expressions  are  numerically  evaluated  for  a  number  of 
points  of  R  determined  by  the  user.  IMSL  subroutines  along  with  custom  written 
subroutines  are  then  used  to  evaluate  the  integrals  in  the  Green’s  functions.  The  non- 
asymptotic  expressions  for  G  a22(R)  and  Gq22(R)  are  evaluated  for  a  predefined  number  of 
points.  The  functions  G^3(R)  and  Gq23(R),  and  the  component  functions  Ib2i(R),  Ib2i(R) 
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and  Ib2i(R)  of  Ge2i(R.C),  and  Ge2i(R.C)  are  then  evaluated.  The  values  calculated  along 
with  their  relative  errors  are  written  to  external  files. 

F.2  SUBROUTINES 

The  following  subroutines  are  called  by  GREENMAIN  or  other  subroutines: 
dzbren  -  IMSL  subroutine  used  to  find  the  root  of  D^X.). 
dqdagp  -  EMSL  numerical  integration  subroutine, 
dqdags  -  IMSL  numerical  integration  subroutine. 

asylimit  -  used  to  find  the  critical  value  for  the  asymptotic  expressions  of 
Ga22(R)  and  Gq22(R).  Inputs:  external  function,  initial  guess  for  stopping 
tolerance  for  result.  Output:  X^. 

chunk  -  Estimates  the  integral  of  a  function  that  oscillates  with  a  Bessel  function 
of  the  first  kind  from  a  to  b  by  integrating  over  finite  intervals  between  a 
and  b  and  summing  the  results.  Inputs:  external  function,  a,  b,  constant  mul¬ 
tiplied  to  variable  of  integration  in  the  Bessel  function,  maximum  relative 
error  allowed.  Outputs:  estimate  of  integral,  estimate  of  absolute  error.  Calls 
subroutine  dqdags. 

bavg  -  Estimates  the  integral  of  an  algebraically  or  exponentially  decaying 
function  that  oscillates  with  a  Bessel  function  of  the  first  kind.  Numerically 
integrates  function  over  an  interval  with  a  user  defined  starting  point  to 
infinity  using  the  Method  of  Averages.  Inputs:  external  function,  order  of 
Bessel  function,  rate  of  decay  for  algebraically  decaying  function,  scale  factor 
if  exponentially  decaying  function,  starting  point,  maximum  relative  error 
allowed.  Outputs:  estimate  of  integral,  estimate  of  absolute  error.  Calls 
subroutine  dqdag  (IMSL). 
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singular  -  Estimates  the  integral  of  a  function  over  the  user  defined  interval  a  to 
b,  with  an  integrable  singularity  singularity  within  the  interval.  Inputs: 
external  function,  a,  b,  free-space  wave  number  k^,  location  of  the 
singularity,  symmetrical  interval  around  the  singularity  to  use  in  cauchy, 
maximum  relative  error  allowed.  Outputs:  estimate  of  integral,  estimate  of 
absolute  error.  Calls  subroutines  dqdagp,  dqdag,  and  cauchy. 
cauchy  -  Estimates  the  Cauchy  principle  value  integral  of  a  function  over  a 
symmetrical  interval  about  the  singularity.  Inputs:  external  function,  location 
of  the  singularity,  distance  around  the  singularity  used  to  define  the 
symmetrical  interval,  maximum  relative  error  allowed.  Outputs:  estimate  of 
integral,  estimate  of  absolute  error.  Calls  subroutine  dqdag. 
infintegral  -  Estimates  the  integral  of  an  exponentially  decaying  function  from  a 
user  defmed  starting  point  to  infinity  by  integrating  over  finite  intervals  and 
summing  the  results  until  the  desired  relative  accuracy  is  reached.  Inputs: 
external  function,  starting  point,  maximum  relative  error  allowed.  Outputs: 
estimate  of  integral,  estimate  of  absolute  error.  Call  subroutine  dqdags. 

F.3  EXTFUNCS 

EXTFUNCS  is  a  file  containing  the  Fortran  implementation  of  the  integrands  for 
the  various  Green’s  functions  and  all  other  non-IMGL  functions  used  in  GREENMAIN 
and  in  other  functions  in  this  file.  The  text  names  of  the  functions  were  given  to  corre¬ 
spond  as  closely  as  possible  to  the  mathematical  function  names  given  in  the  previous 
chapters.  The  following  functions  are  contained  in  EXTFUNCS: 

zsinh(x),  zcosh(x),  and  ztanh(x)  -  double  complex  functions  of  the  respective  hy¬ 
perbolic  functions  where  x  is  also  double  complex. 
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ulb(x),  u2b(x),  and  u3b(x)  -  double  complex  representations  of  UjjjCx),  U2jj(x), 
and  u^jjfx),  respectively,  where  x  is  also  double  complex.  These  functions  re¬ 
quire  the  constant  k2jj  and  k^jj. 

quaddist(i,  n,  rmin,  rmax)  -  double  precision  function  that  returns  value  for  i*^ 
point  in  quadratic  distribution  of  n  points  between  rmin  and  rmax. 

deb(x)  -  double  complex  representation  of  d^{x)  where  x  is  double  precision. 
Function  requires  constants  l^b23’  ^^bl3’  bib-  tjj  and  external  functions 
ulb(x),  u2b(x),  u3b(x),  zsinh(x),  zcosh(x)  and  ztanh(x). 

dmb(x)  -  double  complex  representation  of  d^x)  where  x  is  double  precision. 
Function  requires  constants  £512-  ^b23’  ^bl3’  bib-  *b  external  functions 
ulb(x),  u2b(x),  u3b(x),  zsinh(x),  zcosh(x)  and  ztanh(x). 

imagdmb(x)  -  double  precision  representation  of  imaginary  part  of  dm(x)  where  x 
is  double  precision.  Function  requires  constants  £512*  ^b23’  ^bl3’  bib-  ^b 
and  external  functions  ulb(x),  u2b(x),  u3b(x),  zsinh(x),  zcosh(x)  and 

ztanh(x). 

ngba22(x)  -  double  precision  numerator  of  the  integrand  in  G^2  where  x  is  dou¬ 
ble  precision.  Function  requires  constants  ^523-  *b’  external  functions 
u2b(x),  u3b(x),  zsinh(x),  zcosh(x)  and  IMSL  function  dbsjO. 

realgba22(x)  -  double  precision  function  of  the  real  part  of  the  integrand  in  G^2 
where  x  is  double  precision.  Function  requires  external  functions  ngba22(x) 
and  deb(x). 

imaggba22(x)  -  double  precision  function  of  the  imaginary  part  of  the  integrand 
in  G^2  where  x  is  double  precision.  Function  requires  external  functions 
ngba22(x)  and  deb(x). 

asygba22(x)  -  double  precision  function  the  real  part  of  the  integrand  in  G^2 
minus  the  limiting  value,  where  x  is  double  precision.  This  function  is  used  to 
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find  the  critical  value  for  the  asymptotic  expression  of  Ga^2-  Function 
requires  constants  p.^23  external  functions  u2b(x),  u3b(x), 

zsinh(x),  zcosh(x)  and  deb(x). 

ngbq22(x)  -  double  precision  function  of  the  non-singular  pan  of  the  integrand  in 
'Jq22  where  x  is  double  precision.  Function  requires  constants  |ibi2'  I^b23’ 
^^bl3’  %12’  ^b23’  %13’  ^Ib’  *b’  ^  external  functions  ulb(x),  u2b(x), 
u3b(x),  zsinh(x),  zcosh(x),  ztanh(x),  deb(x),  and  IMSL  function  dbsjO. 

realgbq22(x)  -  double  precision  function  of  the  real  part  of  the  integrand  in  Gq22 
where  x  is  double  precision.  Requires  external  functions  ngbq22(x)  and 
dmb(x). 

imaggbq22(x)  -  double  precision  function  of  the  imaginary  part  of  the  integrand 
in  Gq22  where  x  is  double  precision.  Requires  external  functions  ngbq22(x) 
and  dmb(x). 

asygbq22(x)  -  double  precision  function  the  real  part  of  the  integrand  in  Gq22 
minus  the  limiting  value,  where  x  is  double  precision.  This  function  is  used  to 
find  the  critical  value  for  the  asymptotic  expression  of  Gq22.  Function 

requires  constants  ^^b23’  ^^bl3’  %\2'  %23'  %13’  •’lb*  S  external 

functions  ulb(x),  u2b(x),  u3b(x),  zsinh(x),  zcosh(x),  ztanh(x),  deb(x),  and 
IMSL  function  dbsjO. 

ngba23a(x)  -  double  precision  function  of  the  numerator  qf  the  integrand  in  Ga^ 
where  x  is  double  precision.  Function  requires  r,  external  function  u2b(x)  and 
IMSL  function  dbsjO. 

realgba23(x)  -  double  precision  function  of  the  real  part  of  the  integrand  in  G^3 
where  x  is  double  precision.  Requires  external  functions  ngba23a(x)  and 
deb(x). 
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imaggba23(x)  -  double  precision  function  of  the  imaginary  part  of  the  integrand 
in  G^3  where  x  is  double  precision.  Requires  external  functions  ngba23a(x) 
and  deb(x). 

ngbq23(x)  -  double  precision  function  of  the  non-singular  part  of  the  integrand  in 
Gq23  where  x  is  double  precision.  Function  requires  constants  •^b23’ 
^bl3’  ^bl2’  ^bl2’  ^b23*  ^bl3’  ^Ib’  V  ^  external  functions  ulb(x),  u2b(x), 
u3b(x),  zsinh(x),  zcosh(x),  ztanh(x),  deb(x),  and  IMSL  function  dbsjO. 

realgbq23(x)  -  double  precision  function  of  the  real  part  of  the  integrand  in  Gq23 
where  x  is  double  precision.  Requires  external  functions  ngbq23(x)  and 
dmb(x). 

imaggbq23(x)  -  double  precision  function  of  the  imaginary  part  of  the  integrand 
in  Gq23  where  x  is  double  precision.  Requires  external  functions  ngbq23(x) 
and  dmb(x). 

nibr21(x)  -  double  precision  function  of  the  non-singular  part  of  the  integrand  in 
Ib2i  where  x  is  double  precision.  Function  requires  constants  1^512’  l^b23’ 
^^bl3’  %12’  %12’  %23’  ^bl3’  ^Ib’  V  ^  external  functions  ulb(x),  u2b(x), 
u3b(x),  zsinh(x),  zcosh(x),  deb(x),  and  IMSL  function  dbsjO. 

realibr21(x)  -  double  precision  function  of  the  real  part  of  the  integrand  in  Ib2i 
where  x  is  double  precision.  Requires  external  functions  nibr21(x)  and 
dmb(x). 

imagibr21(x)  -  double  precision  function  of  the  imaginary  part  of  the  integrand  in 
Ib2i  where  x  is  double  precision.  Requires  external  functions  nibr21(x)  and 
dmb(x). 

nibs21(x)  -  double  precision  function  of  the  non-singular  part  of  the  integrand  in 
Ib2i  where  x  is  double  precision.  Function  requires  constants  ltjjl2’  ^b23’ 
^bl3’  %12’  %12’  %23’  %13’  '^Ib’  V  ^  external  functions  ulb(x),  u2b(x), 
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u3b(x),  zsinh(x),  zcosh(x),  cieb(x),  and  IMSL  function  dbsjl. 

realibs21(x)  -  double  precision  function  of  the  real  part  of  the  integrand  in  Ib2i 
where  x  is  double  precision.  Requires  external  functions  nibs21(x)  and 
dinb(x). 

imagibs21(x)  -  double  precision  function  of  the  imaginary  part  of  the  integrand  in 
Ib2i  where  x  is  double  precision.  Requires  external  functions  nibs21(x)  and 
dmb(x). 

nibt21(x)  -  double  precision  function  of  the  non- singular  part  of  the  integrand  in 
Ib2i  where  x  is  double  precision.  Function  requires  constants  e[jl2’  ^bl2’ 
%23’  %13’  ^Ib’  ^b’  ^  external  functions  ulb(x),  u2b(x),  u3b(x),  zsinh(x), 
zcosh(x),  and  IMSL  function  dbsjO. 

realibt21(x)  -  double  precision  function  of  the  real  part  of  the  integrand  in  Ib2i 
where  x  is  double  precision.  Requires  external  functions  nibt21(x)  and 
dmb(x). 

imagibt21(x)  -  double  precision  function  of  the  imaginary  part  of  the  integrand  in 
Ib2i  where  x  is  double  precision.  Requires  external  functions  nibt21(x)  and 
dmb(x). 

F.4  Test  Programs 

The  following  programs  were  used  to  test  the  cauchy  and  bavg  subroutines  on 
functions  with  known  solutions.  Outputs  appear  after  each  program  listing. 

F.4.1  Test  Program  1 

double  precision  £, a, b, pole, delta, errabs,errrel, result, 

&  errest, total, totalerror, exact , error 

c 

external  £,dqdags, cauchy 
c 

a-0.0 
b-3.0 
pole-1.0 
delta-0. 1 
errabs-0.0 
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errrel-le-6 
exact-1 .45059 

call  dqdags (f, a, pole-delta, errabs.errrel, result, errest) 

total-result 

totalerror-errest 

call  dqdags (f , pole+delta, b, errabs,errrel, result, errest) 

total-total+result 

totalerror-totalerror+errest 

call  cauchy (f , pole, delta, errrel , result , errest) 
total-total+result 
totalerror-totalerror+errest 
error-dabs (exact-total) 

write (*, 10) total, exact, totalerror, error 
10  formate  Corrputed  -  ' ,  f  10 . 8, 13x,  ’  Exact  -  ',fl0.8// 

4  '  Error  estimate  -  ' , IpelS. 8, 3x, 'Error  -  *,lpel5.8) 

end 

double  precision  function  f (x) 

double  precision  x 

f-dsin (x) / (x-1) 

return 

end 

Conputed  -  1.45058958  Exact  -  1.45059000 

Error  estimate  -  1 . 46218567E-07  Error  -  4.23679505E-07 

FA. 2  Test  Program  2 

integer  order 

double  precision  dbsjO, alpha, beta, start, errrel, result, 

4  errest, exact, error, r 

external  dbsj0,bavg 

comtton  r 
r-1 . OdO 

order-0 
alpha— 0 . 5d0 
beta-0 . 0 
start -0 . 0 
errrel-le-6 

exact-1 . OdO 

call  bavg (dbsjO, order, alpha, beta, start, errrel, result, 

4  errest) 

error-dabs (exact-result) 

write (*, 10) result, exact, errest, error 
10  formate  Confuted  -  '  ,fl0. 8, 13x,  '  Exact  -  ',fl0.8// 
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4  '  Error  estimate  -  ' ,  IpelS. 8, 3x, ‘Error  -  ',lpel5.8) 

end 

Cotiputed  -  1.00000000  Exact  -  1.00000000 

Error  estimate  -  4 . 44089210E-16  Error  -  1.77635684E-15 

F.4.3  Test  Program  3 

integer  order 
c 

double  precision  f, alpha, beta, start, errrel, result, 

4  errest, exact, error, r 
c 

external  f.bavg 
c 

common  r 
r-1 . OdO 
c 

order-0 
alpha-0 
beta-1 . 0 
start-0.0 
errrel-le-6 
c 

exact-l/dsqrt (1 .0+r**2) 
c 

call  bavg  (f ,  order,  alpha, beta,  start, errrel,  result , 

4  errest) 

error-dabs (exact-result) 
c 

write (*,10) result , exact , errest , error 
10  formate  Conputed  -  ' , f 10. 8, 13x, *  Exact  -  •,fl0.8// 

4  '  Error  estimate  -  ' , lpel5.8, 3x, 'Error  -  ',lpel5.8) 

end 

double  precision  function  f (x) 
double  precision  x,r,dbsj0 
external  dbsjO 
common  r 

f-dexp(-x) ♦dbsj0(r*x) 
return 
end 

Conputed  —  .70710678  Exact  — 

Error  estimate  -  .OOOOOOOOE+00  Error  - 


.70710678 

3.33066907E-16 


FJ  Mathematica  Programs 

The  Mathematica  package  vectorPlot.m  was  used  to  generate  the  vector  plots  in 
chapter  FV.  The  package  polyAvg.m  was  used  to  perform  the  polynomial  average 
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interpolation  necessary  to  obtain  the  data  for  the  vector  plots. 

F^.l  Mathematica  Package  vectorPlot.m 

vectorPlot : :usagea"vectorPlot [data, size]  plots  the 
relative  magnitude  and  direction  of  discrete  data 
points  using  arrows  to  represent  the  output.  The 
input  'data'  has  the  form 

{{{xl,yl,xmag,ymag), {x2, yl, xmag, ymag} ,  . . . 

{ {xl, y2, xmag, ymag} , {x2, y2, xmag, ymag} , ...},.,.} 

xmag  and  ymag  represent  the  x  and  y  components  of  the 
vector  at  position  x,  y.  The  input  'size'  has  the 
form 


{dx,dy} 


where  dx  and  dy  are  the  x  and  y  separations  between 
data  points  in  the  x  and  y  directions,  respectively." 

vectorPlot [data_,  size_] :  = 

Block [ {vector, number, vlength, x, xmin, xmax, y, ymin, ymax, 
mag,  angle,  heads^O  .728869,  tbead=0. 3,  a,b,  c,  d,  arrows} , 
vector=Flatten [data, 1] ; 

(*  Find  the  total  number  of  data  points  *) ; 
numbersLength [vector] ; 

(*  Determine  the  maximum  physical  half  length  of 
a  vector  *) ; 
vlength=44in  [size]  /2//N; 

(*  Extract  the  x-coordinates  and  find  the  minimum 
and  maximum  values  *) ; 
xs^ranspose [vector] [[!]]; 
xmin^fin  [x] ; 
xmax^4ax[x] ; 

(*  Extract  the  y-coordinates  and  find  the  minimum 
and  maximum  values  *) ; 
ysTranspose [vector] [ [2] ] ; 
ymin^Min  [y] ; 
ymaxaMax  [y] ; 

(*  Calculate  the  vector  magnitudes  and  normalize 
to  the  largest  value  *) ; 

magsTable[Sqrt  [vector  [  [i,  3]  ]  ^2+vector  [  [i,  4]  ]  '‘2]  //N, 
{!, number}] ; 
mag^nag/Max [mag] ; 

(*  Calculate  the  angular  direction  of  the  vectors  *) ; 
anglesTable[If [ ( (vector [ [i,3] ]asO) 

&& (vector [ [i, 4] ]»=0) ) , 0, 

ArcTan [vector [ [i, 3] ] , vector [ [i, 4] ] ] //N] , 

(i, number}] ; 
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(*  Calculate  the  head  and  tail  (a  and  b)  coordinates 
of  the  arrow  shafts  *) ; 

a=Table [ {x[ [i] ] +vlength*mag [ [i] ] *Cos [angle [ [i] ] ] , 
y [ [i] ] +vlength*mag [ [i] ] *Sin [angle [ [i] ] ] } , 

{ i ,  niunber }  ]  ; 

b=sTable  [{x[  [i]  ]  -vlength*mag  [  [i]  ]  *Cos  [angle  [  [i]  ]  ] , 
y  [  [i]  ] -vlength*inag [  [i]  ]  *Sin [angle [  [i]  ]  ]  }, 

{ i ,  n\jinber }  ] ; 

(*  c  and  d  along  with  a  determine  the  coordinates 
of  the  arrow  head  *) ; 
c=Table[ {x[ [i] ] +vlength*mag[ [i] ] *head* 

Cos [angle [ [i] ] +thead] , 
y [ [ i ] ] +vlength*mag [ [ i ] ] *head* 

Sin [angle [ [i] ] +thead] } , 

[i,  niJinber}]  ; 

dsTable [ {x [ [i] ] +vlength*mag [ [i] ] *head* 

Cos [angle [ [i] ] -thead] , 
y  [  [  i  ]  ]  +vlength*niag  [  [  i  ]  ]  *head* 

Sin [angle [ [i] ] -thead] } , 

{ i ,  n\jinber }  ] ; 

(*  arrows  contains  the  locations  and  orientations  of 
all  the  graphics  primitives  necessary  to  draw  all 
the  vectors  *) ; 

arrows=Table[ {Line[{a[ [i] ] ,b[ [i] ] )] , 

Polygon [ {a [[i]],c[[i]],d[(i]] }] }, 

{i, number}] ; 

Show[Graphica[ {Thickness [1/800] , 

Line[{{xmin-size[ [1] ] /2,ymin-size [ [2]]/2}, 
{xmax't-size[  [1]  ]  /2,  ymin-size  [  [2]  ]  /2) , 

{xmax+size[ [1] ] /2, ymax+size [ [2] ] /2} , 

{xmin-size[  [1]  ] /2, ymax-t-size [  [2]  ]/2}, 

{xmin-sizei [1] ] /2, ymin-size [ [2] ] /2} }] ,  arrows} , 
AspectRatio->Automatic] ] 


FJ.2  Mathematica  Package  polyAvg.m 

polyAvg: :usage»"polyAvg[data,x]  returns  the  polynomial 
average  for  x  of  data." 

polyAvg [data_, x  ] :» 

Block [{1, il, i?»l,xl, f,outl,out2}, 

IsLen^h  [data] ; 

While [ ( (x>»data [ [i2, 1] ] ) 66  (i2<l) ) , ++i2] ; 
il=i2-l; 

Which  [il=l, 

f-Fit [data[ [Range [1,3] , Range [1,2] ] ] , 

{l,xl,xl^2},xl]  ; 

outl=f/ .xl->x; 

out2soutl, 

i2=»l, 

fsFit [data [ [Range [1-2, 1] , Range [1, 2] ] ] , 
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{1, xl,xl^2},xl] ; 

out2=f/ .xl->x; 
outl=out2. 

True, 

f=Fit [data[ [Range [il-1, i2] ,  Range [1, 2] ] ] , 

{l,xl,xl^2},xl]  ; 

outl=f/ .xl->x; 

f=Fit [data[ [Range [il, i2+l] ,  Range [1, 2]  ]  ] , 

{l,xl,xl^2},xl]  ; 

out2=f/ .xl->x 

]; 

Return [ (outl+out2) /2] 

] 
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>  7^5  7.':  Abstract 

A  theoretical  model  for  the  analysis  of  an  aperture  fed  stacked-patch  microstrip 
antenna  is  presented.  The  mixed  potential  integral  equation  approach  (MPBE)  is  used. 
The  aperture  is  closed  by  using  opposing  magnetic  currents  on  each  side  of  the  ground 
plane.  The  necessary  Green’s  functions  associated  with  the  vector  and  scalar  potentials 
are  evaluated  in  the  spatial  domain  using  stratified  media  theory.  The  Green’s  functions 
are  expressed  as  Sommerfeld  integrals.  A  method  of  moments  technique  to  solve  for  the 
currents  of  the  antenna  is  outlined.  Basis  and  test  functions  are  found  to  use  with  the 
Green’s  functions  in  the  MPIEs  to  form  a  solution  matrix.  No  actual  solutions  for  the 
currents  are  calculated. 


The  Sommerfeld  integrals  in  the  Green’s  functions  are  analyzed  to  determine 
their  characteristics.  These  characteristics  include  complex,  oscillatory  integrands;  sin¬ 
gularities;  surface  waves;  and  semi-infinite  integration  intervals.  Several  numerical  inte¬ 
gration  techniques  to  deal  with  these  characteristics  are  developed.  Both  Fortran  code 
and  Mathematica  packages  written  to  implement  these  techniques  are  discussed. 

Example  calculations  for  several  of  the  Green’s  functions  are  accomplished.  The 
cut-off  frequencies  for  the  surface  waves  are  evaluated  and  it  is  shown  that  with  a  proper 
choice  of  material  parameters  only  one  surface  wave  mode  will  propagate.  The  Green’s 
functions  are  then  evaluated  accurately  and  efficiently  with  sample  results  provided. 


Ideas  for  continued  research  and  new  applications  are  discussed. 
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